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Abstract. In the first part of this paper wc construct an algo- 
rithm for implementing the discrete wavelet transform by means of 
matrices in SO2 (R) for orthonormal compactly supported wavelets 
and matrices in 5Lm(R), > 2, for compactly supported biorthog- 
onal wavelets. We show that in 1 dimension the total opera- 
tion count using this algorithm can be reduced to about 50% of 
the conventional convolution and downsampling by 2-operation for 
both orthonormal and biorthogonal filters. In the special case of 
biorthogonal symmetric odd-odd filters, we show an implementa- 
tion yielding a total operation count of about 38% of the conven- 
tional method. In 2 dimensions we show an implementation of this 
algorithm yielding a reduction in the total operation count of about 
70% when the filters are orthonormal, a reduction of about 62% 
for general biorthogonal filters, and a reduction of about 70% if 
the filters are symmetric odd-odd length filters. We further extend 
these results to 3 dimensions. 

In the second part of the paper we show how the S'02(R-)- 
method for implementing the discrete wavelet transform may be 
exploited to compute short FIR filters, and we construct edge map- 
pings where we try to improve upon the preservation of regularity 
due to conventional methods. 

In the third part of the paper we consider the problem of 
discriminating two classes of radar signals generated from some 
number of point sources distributed randomly in a bounded plane 
domain. A statistical space-frequency analysis is performed on 
a set of training signals using the LDB-algorithm of N.Saito and 
R.Coifman. In this analysis we consider several dictionaries of 
orthonormal bases. The resulting most discriminating basis func- 
tions are used to construct classifiers. The success of different 
dictionaries is measured by computing the misclassification rates 
of the classifiers on a set of test signals. 
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CHAPTER 1 



Preliminaries and Concepts 

1. A library of bases for L'^(R) 

Let L'^(R) denote the space of all square integrable functions of a 
single real variable. We present some well-known bases for this space. 

1.1. Wavelet bases. Expressed shortly, these bases consist exclu- 
sively of all dyadic dilations and integer translations of a single mother 
wavelet ip with zero integral, that is 

/ G L\R) ^f=Y^ d^. , where i^j^x) = 2^/2^(2^x - k), 

with the series converging in sense. In this paper we will only be 
concerned with real and compactly supported wavelets. In practice, ip 
will possess some localization in both time and frequency, and the basis 
will have to be stable, thus asserting the existence of positive numbers 
A and B such that 

and ip will possess some vanishing moments, that is 

J i^{x)x^dx = 0, < Z < M, M > 1. 

We can further separate this family of bases into two subfamilies: 

• Orthonormal Wavelet Bases. The mother wavelet is orthonor- 
mal to all its dyadic dilations and integer translates, thus assert- 
ing 

• Biorthogonal Wavelet Bases. The mother wavelet is biorthogo- 
nal to all dyadic dilations and integer translates of its dual ■0, 
yielding 

1pj,k,'lpj',k'^ = Sj,fSk,k' , dj^k = (^f,i^j,k^ ■ 
9 
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The ipj^k are called the analysis wavelets, while the ipj^k are called 

the synthesis wavelets. 
Wavelet bases are intimately connected to the concept of a "Multires- 
olution Analysis" (MRA). 

Definition 1.1. A Multiresolution Analysis is a nested sequence 
{Vj}j(zz of closed subspaces o/L^(R) satisfying 

1. V, c V,+ i ,jez. 

4. f{x) e Vj =^ f{x -k) eVj ,keZ. 

5. fix) G V, =^ f{2x) G 

6. There exists a "scaling function" G Vq such that {0j,A;}fc6Z con- 
stitute a Riesz basis (= stable basis) for Vj , j ^ Z. 

The great triumph of MRA is the following result which proof can be 
found in [jl|: 

Theorem 1.1. Given a MRA, there exists a wavelet ijj e ViD 
such that Spaiij j^ipj^k = L'^(R.)- The wavelet can be chosen to have 
compact support. 

Remark. The wavelets possessing the maximum number of vanishing 
moments compatible with their support width, are called Daubechies' 
wavelets. It is shown in that a wavelet with vanishing moments 
has support width 2N — 1. It is also possible to assign vanishing mo- 
ments to the scaling function 0, (except a zeroth vanishing moment). 
Wavelets where both the wavelet and the corresponding scaling func- 
tion have vanishing moments are called coiflets. 

From the first MRA property we get the "scaling iteration equation" 

= 2^/2 5^ if (/ - 2A;)0,+i,z , H{1) = {(t>Ai,i) ■ (1-1) 

If the {0o,fe}fcgz form an orthonormal basis for Vq (if not we can carry 
out an "orthonormalization trick", see for details) we get an or- 
thonormal wavelet basis from the MRA by defining the wavelet space 
Wj at scale j as the orthogonal complement of Vj in V^+i, that is 

Vj+i = Wj ®V.j. (1.2) 

Now, the properties of the scaling spaces Vj together with the definition 
(|1.2|) yields the orthogonal decomposition 

L2(R) = 0iy„ (1.3) 
and the "wavelet scaling equation" 
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V-.^fc = 2^/2 - , Gil) = iij, 01,^) . (1.4) 

The vanishing integral of ip imphes a non-vanishing integral for 0. 
To see this, assume / = and consider the characteristic function 
X[-K,K] G L^- By choosing K large enough, we get a contradiction to 
the fact that {4>j^k}j,kez constitute a Riesz basis for L^. 

Since En H(n) = / ^ and En G{n) = J i/j = 0, the H{n) 
make up a lowpass filter, convolving it with x G R" produces weighted 
local averages, while the Qn make up a highpass filter, convolving it 
with X G R" produces details by catching up local oscillations. Filters 
with finite support width (i.e finite number of nonzero coefficients) are 
called "finite impulse response filters" or simply FIR filters. 

By multiplying equation (|1 . 1| ) on both sides by (pj^k and integrating 
the result, using orthonormality of to its integer translates, we get 
orthonormality of the lowpass filter to its even translates 



5o,fc = J]^(n)iJ(n + 2A;). (1.5) 

nGZ 

In |I| it is shown that the highpass coefficients G{n) are uniquely de- 
termined by the lowpass coefficients H{n) up to 



G{n) = {-iy'H{2N +l-n) , modulo phase and N e Z. (1.6) 

The finite support of the wavelet ip and the filters {H{k)}kez, {G{k)}kez 
are thus seen to be simple consequences of the finite support of the scal- 
ing function. We can also see that by relation ( p,.6| ) the lowpass and 
highpass filters have the same length, moreover, it has to be even by 
relation (|1.5|). 

We remark that by Fourier transforming equation (11 . 11) and iterat- 
ing the result, we get a (infinite) product formula for 0, given by 



oo 

— n^o(2-^o 



L-l 



V2 



n=0 



thus (and consequently ip) is uniquely determined by the filter coef- 
ficients H{n). 

Now, given the projection Pjf of a function / on some scaling 
space Vj, there is an easy way of computing the projection of / on all 
"coarser" spaces {Vj, Wj}j^j without integrating. From ( |1.1] ) and ( |1.4| ) 
we obtain a set of recursion relations: 
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(/, = 2^/2 J2 H{1 - 2k) (/, 0,,,) , (1.7) 

if, = 2^/2 ^ G(/ - 2k) if, 0,,,) . (1.8) 

We observe that at each scale j the projection Pjf is decomposed into 
two parts: 

• Local averages. The lowpass coefficients = (/, (pj,k) are the 
resuh of passing {c^~ }/gz through the "lowpass channel" deter- 
mined by (|1.7|) . This operation can be described as "convolution 
and downsampling by 2": First convolve with {H{—k)}k(^z, then 
throw away the odd numbered indicies. This operation is often 
called "lowpass filtering". With the convention c = {c„}„gz, we 
can write a (infinite) matrix equation 

c^-i = He'', U{iJ) = H{i - 2j), < i,j <oo. (1.9) 

• Local oscillations. Similarly, the highpass coefficients dl. result 
from passing {c^^^j^ez through the highpass channel determined 
by ( |1.8| ) and is called "highpass ffitering" . We write 

rfi-i = Gc^' , G(i, j) = G{i - 2j), 0<ij <oo. (1.10) 

The decomposition of Pjf into averages and details at different 
scales is illustrated in Figure 0. This decomposition is the discrete 
wavelet transform in dimension 1. 




Figure 1. The "logarithmic tree" of a wavelet decom- 
position in M levels. 



The inverse operation is in this orthogonal case given by the transpose 
matrices: 

4 = Hy-\k) + G'd^'\k) 

= 2^'^^H{k-2l)(^~^ + 2^/^^G{k-2l)d^-\ (1.11) 
iez «ez 

This operation is illustrated by simply reversing the arrows in Figure 
|1| and replacing H, G by their transposes. 
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We emphasize that the expansion of a function into a wavelet basis 



is fast: From the equations (|1.7| ) and (|1.8| ) we get that a complete 
wavelet analysis of a sequence of length N has a total operation cost 
of no more than 4L ■ N multiply-adds, where L is the support width 



of the filters. Equation (|1.11|) yields the same bound on the operation 
cost for reconstruction. 

In summary the FIR filters {H{n)}nez, {G{n)}nez make up a two- 
channel orthogonal filter bank with perfect reconstruction. 



If is not orthogonal to its translates, the sums in ( |1.2| ) and (jL 
are direct rather than orthogonal and the wavelet basis is no longer 
its own dual. The construction of a biorthogonal Riesz wavelet basis 
is shown in Thus we get two MRA^s: The analysis MRA that is 
built up from the analysis scaling functions (pj^k and the corresponding 
wavelets ipj^k, and the synthesis MRA built from the synthesis scaling 

functions 4>j^k and their wavelets ipj,k- 

The relations between the basis functions and filters from the two 
MRA^s are given below. A simple picture of the biorthogonality rela- 
tions is shown in Figure ^ 



Filter relations: 



;i.l2) 



G{n) = {-IY+^H{2N + 1 - n), 

G{n) = (-l)"+^if(2Ar + l-n),N eZ. 



:i.l3) 



Biorthogonality relations: 



'ipj,k,'ipj',k'j = Sj,j'Sk,k' {(po,k,(pQ,k'j = Sk,k'- (1-14) 

The scaling equations for the biorthogonal scaling- and wavelet func- 
tions and the expansion and reconstruction equations are immediate 
generalizations of the corresponding equations in the orthogonal case. 
We note that the lengths of the filters H and H need not be equal. 

When working in L^(R"), that is with square integrable functions 
of n variables, it is possible to generate a wavelet basis for this space by 
taking tensor products of n (possibly different) MRA's in L^(R). This 
results in 2" — 1 different mother functions, the dilations and translates 
of which make up what we call a tensor wavelet basis in dimension n. 
The iteration is done on the pure lowpass band, that is the coefficients 
originating from the pure tensor scaling function 0(xi) ® 0(2:2) ® ■ ■ ■ ® 

(j){Xn). 
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Figure 2. The scaling- and wavelet spaces Vj, Wj and 



their duals Vj, Wj. 



1.2. Wavelet packet bases. These bases are particular linear 
combinations of elements in a wavelet basis. It follows that the ele- 
ments in a wavelet packet basis inherit some of the properties of the 
wavelets they are made of, such as compact support and some local- 
ization in both time and frequency. We will see that the wavelet basis 
is included special case of wavelet packet bases. 

Following p|, we briefly describe the construction of these bases. 
Define 



iez 

i,s,fAt) = 2-^/V/(2-^t -p); Af = Span7/;o,/,p, 



p 



ax 



(t) = 2-'^/^x{t/2); aAf = {ax : x E A/}, (1.15) 



where if, G is a pair of conjugate FIR filters and 'ip.s,f,p is called a wavelet 
packet of scale index s, frequency index / and position index p. We 
immediately have that a'^Af is the closure of Spanpipsj,p- Exploiting 
the natural one-to-one correspondence 



'sj ' ' (^"Af- Is J 



s s + l 



2/' 2/ 

we have the following result for which proof we refer to |2|: 

Theorem 1.2. If I is a dyadic cover o/R"*", then 
{IJ^ jCr^Aj : Isj G X} = i^^(R) and if I is disjoint, the wavelet packets 
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{'4^s,f,p '■ Isj ^l-iV ^ Z} jorm a basis for L^(R). Furthermore, if the 
filters H , G are orthogonal, this basis is an orthonormal basis. 

To keep the same filtering formulas for sequences as for functions, make 
the definitions: 



Then it is easy to verify the following recursion relations: 



s+l,2f,p 
s+l,2/+l,p 



HX 



G\ 



:i.i6) 



Thus, identifying cr'^Af with fisj, the expansion of a function into a 
collection of wavelet packet bases can be illustrated as in Figure 0. We 
note that every disjoint cover of the top box in this figure by smaller 
boxes from levels below, yields a new basis for the space floo- 



Scale 





^1,0 






^\ 


H / 




^2,0 


^2,1 


^2,2 


^2,3 



G H \G H / \G Hf \G 



^3,0 


^^3,1 


^^3,2 


f^3,3 


^3,4 


^3,5 


^3,6 


^^3,7 



Frequency 

Figure 3. The binary tree of a wavelet packet decom- 
position in 4 levels. 



By comparing Figure |I] and Figure |^ it is easy to see that the logarith- 
mic tree of expansion into a wavelet basis is a subtree of the binary 
tree of expansion into a collection of wavelet packet bases. We conclude 
that the wavelet basis is included in the collection of wavelet packet 
bases. 

For reconstruction one needs the dual basis of the particular wavelet 
packet basis. It is shown in that the wavelet packets defined by the 
dual filters H' , G' are the duals of the wavelet packets defined by H 
and G. 
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Given a sequence of length N, a rather crude estimate given in 
0] shows that a wavelet packet analysis of the sequence will provide 
more than 2^ different bases. The recursion relations ( |1.16| ) makes the 
expansion into a wavelet packet analysis cheap: The total operation 
cost will be no more than L ■ N log2 N multiply-adds, where L is the 
support width of the filters. 

1.3. Smooth local trigonometric bases. The huge drawback of 
the Fourier basis for applications in signal analysis is its "non-locality" 
in time, thus it is impossible to relate specific frequencies of a signal to 
specific time or space windows using this basis. Also, this basis only 
span the space of the periodic L'^(R) functions. What we would like is 
a trigonometric basis for L'^{R.) with the following desirable properties: 

• Orthogonality. 

• Smoothness: The basis elements should possess at least a few 
continuous derivatives. 

• Localization in space: Each basis element should have finite sup- 
port. 

• Efficiency: There should exist a fast algorithm (i.e. N logg A^) 
for expansion into the basis. 

Such bases do exist, we refer to for a thorough exposition on the 
subject. We will briefly outline the construction of a Local Sine Basis 
and prove that this basis enjoys the properties above. Define the bell 
functions {bk{x)}kez with the following properties: 

J2blix) = l, VxeR, 

supp{hk) = Ik = («fc - Cfc, Ofc+i + efc+i), 
bkbk+2 = 0; bk-ibk is even around a^. (1-17) 
A schematic picture of these bell functions is shown in Figure |^. 



1 1 




R 



Oik ttfc+l Olk+2 



Figure 4. The bell functions. 

Define 

s,k{x) = ^'^^^ sin Uj + 1/2) ^""^ ) . (1.18) 
(afc+i-«fc)^/^ V ak+i-akj 
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Then we have the following result: 
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Theorem 1.3. The sj^k form a smooth orthogonal basis for L^(R). 
Furthermore, this basis enjoys all the desirable properties listed above. 

Proof: The smoothness and localization properties of this basis are 
immediately clear. To prove orthogonality we only need to consider 
{sj^k, Sj',k') when \k — k'\ < 1 because of the support property of the 
bell functions. Define 



~/\ 1 f / , X — ak \ 

(afe+i - afejV^ V oik+i - Oik J 

For k — k' — 1 we get 



{sj^k,Sj',k~i) ^ / bk-i{x)bk{x)sf^k-i{x)sj^k{x)dx ^0. 

To see this, observe that since bk~ibk is even around ak, the integrand 
is an odd function around the center of integration, thus the integral 
equals zero. For k — k' we get 



/■afc+i+^fc+i 

{sj,k,Sf,k) = / bl{x)sy,k{x)sj,k{x)dx 

Sji^k{x)sj^k{x)dx = 6jj'. 



To see that the first integral equals the second, observe that Sf^k'sj^k is 
even around ak and b\{x) = 1 — &^,(afe + — x), x E[ak — ek,Oik + €k\. 
The second integral is an elementary calculation. This shows that the 
Sj^k are orthonormal over all frequencies j and translates k. 

Now, let / e -^^(R) and consider (/, Wc split the integral into 
three parts and change variables on the leftmost and rightmost integrals 
to obtain a formula for the inner product involving integration only on 
the interval We may express this operation by saying that 

the part of / that lives in [ak — e^, 0;^] U [aj^+i, ak + Cfe+i] is folded into 
the interval [ccfe, CKjk+i]. 
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{f,Sj,k) = / f{x)hk{x)sj^k{x)dx 

J oik-ek 

(■)+/ (■)+/ (■) 

fOik + l 

bk{x)f{x)sj^k{x)dx 

' 0!k 

bk{2ak - x)f{2ak - x)sj^k{x)dx 

r^k+i 

+ / fefc(2afc+i - x)f{2ak+i - x)sj^k{x)dx 

J oik 

"O^k + l 

f\x)sj^k{x)dx, 

'oik 

fix) = bk{x)f{x)-bk{2ak-x)f{2ak-x) 
+ bk{2ak+i- x)f{2ak+i- x). 

To implement the local sine basis, we place the gridpoints at half- 
integers, and evaluate the integral by DST — IV, meaning a discrete 
sine transform of type IV. This transform enjoys the nice property of 
being its own inverse, and has a N logg implementation. We refer to 
0] for a proof of these facts. 

It only remains to prove completeness. We have (/, Sj^k) = {fK^j,k)- 
Furthermore, Span^Sj^fc = L'^{\oik,Oik+i)) because of the completeness 
of the Fourier basis in L^([0,27r)). This yields 



oo 

bk{x)f{x) = ^ (/^Sj-fc) &fc(a;)sj,fe(x) 

j=0 

oo 

j=0 

Pk{f){x) = bl{x)f{x) 

- bk{x)bk{2ak- x)f{2ak- x) 

+ bk{x)bk{2ak+i- x)f{2ak+i- x). 

Summing Pk{f) over all G Z, it is easy to see that all the terms 
cancel except the b\{x)f{x) terms. Thus we get 
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feez j=o feez Vfcez / 

this confirms the completeness of the Sj^k- ^ 

We note that by replacing sines by cosines and DST—IV by DCT— 
IV in the construction above, we get a local cosine basis for L^(R). 

1.4. Local sine/cosine packet analysis. It is immediately clear 
how the notion of a wavelet packet analysis can be transfered to this 
case of local trigonometric bases. We restrict ourselves to bells of a 
fixed length at each level of decomposition. Thus, we replace each bell 
from the level above by two child bells of half the length of the parent 
bell. A schematic picture of this analysis is shown in Figure ^. 




Figure 5. Dyadic local sine/ cosine analysis in three levels. 

It is easy to see that in analog to the case of wavelet packets, every cover 
of the original interval by our special bell functions corresponds to a 
(unique) basis for this interval. To get back to the original coordinates, 
simply perform the DST — IV or DCT — IV (they are their own 
inverses) on each subinterval, then "unfold" the result back into the 
folding regions [au - ek,ak + e^]. 

For a input of length N, the number of levels in a Local Sine/ Cosine 
Analysis cannot exceed log2 A^, thus the total operation cost is bounded 
by A^(log2 A^)^ multiply-adds. 

2. Signal classification using local feature extraction 

The notation and ideas presented below are due to [^, to which 
we refer for a more thorough review of this subject. We say that a 
wavelet packet analysis or a local (co)sine packet analysis constitute a 
dictionary of bases for L^. We call a collection of dictionaries a library 
of bases for L^. 
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2.1. Local feature extraction. We define tlie problem of signal 
classification as the construction of the map c 

c : 5 C I — >C = {1,2,...,C}, 

where S is the set of all signals under consideration and C is the set 
of all relevant class names. We call S the signal space, C the response 
space and c a classifier. Since the dimension n of the signal space is 
normally very large compared to the dimension of the response space, 
it is important to extract only the relevant features of the signals to 
obtain an efficient and accurate classification. In other words, we want 
to map the information relevant to our problem into a few coordinates 
and ignore all the rest. To achieve this, define the map / 

/ : 5 I — ^ C R^ k<n, 

T is called the feature space and / the feature extractor. Thus, once 
we have /, we only have to construct c : i-^ C, which in general is 
much easier because k « n. 

To construct the feature extractor we use a training data set r = 
{^i)Ci}f=i C S X C of N training signals Sj and their class names q. 
We will denote by Ni the number of training signals belonging to class 
i, so that = A''i + ■ ■ ■ + Nq. Once / and c have been constructed, 
we measure the mis classification rate test data set r' = 

{s^, c'f}fLi which has not been used in the construction of / and c, by 

1 

Terror = J^^^i^'i^ ^(s-)), 5(0) =0, 5{x) = 1, 0, 
i=l 

We will use feature extractors on the form / = O*^'^-' o \[' where ^ G 
S02{n) with column vectors ^J^k,i "with the correspondence 

SpanwJ'fc, = fij-fc, 
I 

J = 0, J, A; = 0, 2^' - 1, / = 0, 2" - 1, (1.19) 

where flj^k are the subspaces of a wavelet packet or a local sine/cosine 
analysis. G'-'^^ : S i — > is a selection rule picking out the k most 
relevant coordinates from n coordinates. 

Figure |] illustrates two examples of C = {1,2} using some \E' and 
some Q^'^\ The classifiers c are indicated in the figure. 

There are several possible choices for a classifier c, such as Linear 
Discriminant Analysis (LDA) and Classification and Regression Trees 
(CART), we refer to for a review of these, as we will not make use 
of them in this paper. 
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Figure 6. Example of sample plots of the two most 
discriminating coordinates in the two-class case. 

2.2. Entropy and discriminant measures. Given a norm || • 
we define the entropy Hr of a sequence x as 

I ir II'' 

^r(x) ^-Y^ -p^ log2 -p^, 1 < r < oo. 

■X W r -X r 

Hri^) measures the degree of order in the sequence x, that is the 
information cost in describing x. 

In the two-class case, we need a discriminant measure d that mea- 
sures how differently two sequences are distributed, in other words the 
relative entropy of the two sequences. In this paper we will only use 
d = m„ defined as 



mp(x, y) = ||x - y||^ = ^(x, - y,f , (1.20) 

i 

with p = 2 in most cases. 

In the general case of C classes, we define the discriminant measure 
of C sequences as 



d({x('^)}f=,) ^ Xi E di^^^^^'^)- (1-21) 

i=l j=i+l 

We see that once we have a discriminant measure, we are capable of 
evaluating the power of discrimination of the different subspaces in any 
dictionary V. This indicates that it is possible to choose a basis in 
for the set <S of signals with the property that no other basis in the 
dictionary will discriminate more between classes. Furthermore, the 
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discriminant measure should be additive to ensure a fast computational 
algorithm. 

Definition 1.2. A discrimant measure d is said to he additive if 

2.3. The local discriminant basis algorithm. Given a discrim- 
inant measure d, how do we best evaluate the power of discrimination 
in each subspace of a dictionary V, and how do we select the most 
discriminating basis? The result should ideally only depend on the 
characteristic features of each class. The "local discriminant basis al- 
gorithm" (LDB-algorithm) developed in |Q yields a particular basis 
called a local discriminant basis, LDB for short, as described below. 

Definition 1.3. Let be a set of training signals belonging 

to class c. Then the time-frequency energy map of class c, Tc , is a table 
of real values specified by {j, k, I) as 

Uj,kj) = J2{<.y4^) (1-22) 

1=1 i=l 

forj = 0,...,J, A; = 0,..., 2^-1, / = 0, 2--^- - 1. 
For notational convenience, define 



2"-J_l 

d({r,(j,fc,-)}f=i)= Yl rf(riUfc,/),...,rc(j,fc,0), 

Bj^k = Span ^j,k,i C flj^k, 

0<«<2"-^-l 

Aj,k = LDB|^^.^ , 

^j,k = discriminant measure of Qj^k- 
Then we have the following algorithm: 

Algorithm 1.1. (The Local Discriminant Basis Selection Algo- 
rithm). Given a training dataset r consisting of C classes of signals 

{{4^}r=i}^=i, 

Step 0: Choose a dictionary T> of orthonormal wavelet packets or local 
sine/cosine packets and specify the maximum depth J of decomposition 
and an additive discriminant measure d. 

Step 1: Construct time-frequency energy maps for c = 1,...,C. 

Step 2: Set Aj^k = Bj^k and Aj^k = ci({r,( J, A;, ■)};?=i) for k = 
0,..., 2-^-1. 
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Step 3: Determine the best subspace Aj^k for j = J — 1,...,0, k = 
0, 2-^ — 1 by the following rule: 

Set A,-, = f/({r,(J,A;,-)}f=i). 

If ^j,k > '^j+l,2k 

+ A 

then Aj^k = Bj^k, 

else Aj^k = ^i+i,2fc © Aj+i^2k+i and set Aj^k = ^j+i,2k + Aj+i^2fc+i- 

Step 4: Order the basis functions by their power of discrimination 
(explained below). 



Step 5: Use k < n most discriminating basis functions for constructing 
classifiers. 

We note that Step 3 is fast, 0{n), since d is additive. In Step 4, 
we evaluate the power of discrimination of a single basis function "Wj^k,i 
by (i({rc(j, fc, That the basis obtained by the LDB-algorithm 



indeed has the desired property, is stated in Proposition for which 
(simple) proof we refer to 

Proposition 1.1. The basis obtained by the LDB-algorithm maxi- 
mizes the additive discriminant measure d on the time-frequency energy 
maps {Pc}^! among all the bases in the dictionary T>. 
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CHAPTER 2 



Factorization of the Discrete Wavelet Transform in 
Dimension 1 Using Elements in S02(Sl), SLm(R) 

1. Introduction 

1.1. Notation. If T is a linear transformation, tlien T^,T~^ will 
be the transpose of T,T~^, respectively. If x e R", then S{k)x{j) = 
x(j - k) and ax(j) = x(-j). If a : Z" ^ R", 6 : Z" ^ R", then 
a * b{k) — X]jGZ^O)^(^ ~ j)- ^"^^ convenience we will occasionally 
denote Daubechics' shortest filters of length L, corresponding to an 
orthonormal compactly supported wavelet possessing a number of L/2 
vanishing moments, by Hf for lowpass and for highpass. We will 
denote the operation of "convolution followed by downsampling by 2" 
by *2- Using similar notation, in the biorthogonal case we will denote 
the analysis lowpass- and highpass filters by H'l, respectively, and 

the synthesis lowpass and highpass filters by G\, where L and L 
are the lengths of the lowpass analysis- and lowpass synthesis filter, 
respectively. 



1.2. The idea. We can present the main idea in our construction 
in just a few words. Consider a conjugate pair of FIR analysis filters 
H, G. The filtering of a vector in R", n even, into lowpass and highpass 
coefficients using this pair of filters for the operation *2, can be arranged 
as a mapping induced by a linear operator Qn{H, G) : R" ^ R" defined 

by 



e„(//, G) : X e R" ^ (cj, dl cl, d}, • • • , 4/2-1, </2-i) ■ (2-1) 

Given a vector x G R", we consider the following linear operators 
from R" to R", defined by 



F2(a) : X ^ |M2(a) 
S02{K) 3 M2{a) 



X2j 
1 



VTT 



a" 



n/2-l 

1 —a 
a 1 



(2.2) 
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, ^ ^ n/m—\ 



,m>2, 



j=0 



1 if i — j 

SL^(R)3 MMihj)- { -« if^ = l,i = m (2.3) 

otherwise. 



SLs(R) 3 Ms{a) = 1 , (2.4) 




5 



5L2(R) 9 M2(a) = ( i "1, (2.5) 



l + aP\P 1 

Remark. In case the number n/m,m > 2 is non-integer, just extend 
X e R" by periodizing or mirroring or by zeros. 

We will show how the operator ©n(i?, G) can be factored into: 

• A composition of operators -p2(«j) and S{kj), j in some finite 
index interval, in the case where H,G is a. conjugate orthogonal 
pair of filters. 

• A composition of operators F^^{aj), F^{aj) and S{kj), with m, 
j in some finite index intervals, in the case where if, G is a pair 
of biorthogonal conjugate analysis filters of odd-odd lengths or 
odd-even/even-odd lengths. 

• A composition of operators F^^{aj), F^^{aj) and S{kj), with j 
in some finite index interval, in the case where if, G is a pair 
of biorthogonal conjugate symmetric analysis filters of odd-odd 
lenghts. 

• A composition of operators F~^{aj), F^^(aj), F2'^(aj), F2 *(<^j) 
and S{kj), with m, j in some finite index intervals, in the case 
where H,G is a conjugate pair of biorthogonal filters of even-even 
lengths. 

We will start considering Daubechies/Coifiet filters, and then ex- 
tend our results to biorthogonal filters. 
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2. Factorization of the orthonormal wavelet transform. 

Let H'l, Gj^ be an orthonormal pair of FIR filters. A linear orthog- 
onal map on is some matrix M(a), a G R in the one-parameter 
family (R). If we can represent 0„(if^, ) by some composition of 
the orthogonal maps F2{aj) as defined in (|2.2|), we have an alternative 
way of computing the discrete wavelet transform in the orthonormal 
case. 

Defining the normalization factor R^^ by 



Rm = X{{l + a])-'^, (2.6) 
i=i 

then because of the linearity of the maps -F(aj), we can avoid repeated 
multiplication by square-root factors, instead normalizing in the last 
step by scalar-multiplication by some Rm- 

From the orthogonality relation (|1.5|) we see that the operation of 
lowpass filtering when applied to the lowpass filter itself, reduces the 
filter to 1 point. Obviously, we have to construct a stepwise procedure 
to reduce the length of the lowpass filter, using some set of matrices 
M2{aj). We define the first step to be the operation 



Ht^F^MHt (2.7) 



where we claim 



Hl{0) ^ 0, which implies ai = — (2.8) 



Then by ( |1.5D we have 



= HtiO) . Ht{L - 2) + Htil) . HtiL - 1) 
^ HtiO) HtjL - 1) 
Hlil) Hi{L-2) 
=^ Hi{L - 1) ^ under ^2(01). (2.9) 

We see that the number of coefficients in the filter F2{ai)H'[ is 2 less 
than in Hf. It is now clear how to proceed: To simplify notation we 
write 



W2{ak) = F2{ak)S{-{k - 1)). (2.10) 
W2(-) is an orthogonal operator. Then we define the fc'th step in our 



procedure, 1 < k < L/2 as the operation 
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W2iak-i) ■ ■ ■ W2{ai)Hi ^ l^2(«fc)W^2(afc-i) ■ ■ ■ 1^2(«i)i^L- (2.11) 

where we claim W2{ak) o ■ ■ ■ o PV2(ai)-ff^(0) t— > 0. 

Because every (oj) is a linear orthogonal operator, the orthogo- 



nality relation (|1.5| ) is preserved in each step. Thus each step but the 
last will reduce the length of the filter at the previous step by 2, the 
last step reducing a filter of length 2 to a filter of length 1. Our con- 
struction is illustrated in Figure below for the filter H^. Each cross 
in Figure m represents the mapping induced by M2{a) on the lower pair 
of points resulting in the upper pair of points. 

M2(Q3) 

1. -0 




Figure 1 . Orthogonal mapping of the lowpass filter Hq 
to 1 point. 

Since the highpass filter is the alternating flip of H'l (see relation 
(|1.6|) ), it is easy to see that applying the operation defined in to 

will reduce it to a 1-point-filter. Because of the orthogonality of 
lowpass channel to highpass channel, we see that 



Ht ^ (1,0) ^Gi^ (0,1). 
Thus we get the following result: 

Theorem 2.1. Given a pair of orthonormal FIR filters H'[, G'l of 
length L, the operator Qn{H'[,G'l) defined in 4 j may be decora- 
posed as QniHtGi) = S{1 - L/ 2)W2{a l /2)W2{a l /2-1) ■■ - W^ia^), 
with W2{aj), I < j < L/2 defined as shown above. 

Proof. The theorem was proved by construction. The argumen- 
tation below is just a formalization of this construction. 

cl = {^*2aHf){k) = Y,Ami{j-2k) = {S{2k)Ht^) 

iez 

= (W^2"'(«i) ■ ■ ■ W^\aL/2)W2{aL/2) ■ ■ ■ W2{ai)S{2k)Ht, x) 
= <W^2(«L/2) ■ ■ ■ W2{a,)S{2k)Ht, W^\aL/2) ■ ■ ■ W^\ai)^) 
= {S{j, 2k - {L/2 - 1)), W^2K/2) ■ • • H^2(«i)x> 
= W^2K/2)---W^2(ai)x(2A;-(L/2-l)). (2.12) 
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4 = {^*2CTGi){k) = J2^U)G{j-2k) = {S{2k)Gi,^) 

= {W,-\a,) ■ ■ ■ W,-\aL/2)W2{aL/2) ■ ■ ■ W2{ai)S{2k)G'i, x> 

= {W2{aL/2) ■ ■ ■ W2{a,)S{2k)Gl W,-\aL/2) ■ ■ ■ W,~\a,)^) 
= 2k + 1- (L/2 - 1)), W2iaL/2) ■ ■ ■ W^2(ai)x) 

= W2{aL/2)---W2{aM2k + l-{L/2-l)). (2.13) 

□ 

A schematic picture of the algorithm that results from the factor- 
ization of G„(if^, G'l) is shown in Figure |^ below. 



Cq cIq c\ d\ c\ d\ 



■Rl/2 -»- • • 

W2(a3)— »- 
1^2(a2)— »- 

W2{ai)^ /\ ^ ^ /\ 
x(0) 




/2-1 



x(n - 1) 



Figure 2. Filtering-diagram for Hq,Gq. 



The "tilde" c's and d's at the edges of the filtering-diagram above 
can be obtained by mirroring x about its endpoints or periodizing. 
Counting operations, we see that we have to do 2 multiplications and 
2 additions pr. pair of points for each 14^2 Writing 



1 

Ilr^iHtGi)^ n ^^K), (2.14) 

j=L/2 

then by construction S{1 — L/2)n„ = 0„, but it has a more efficient 
implementation as shown in Table |^. [] 



Operation 


n„ 


e„ 


tl mult's 


(L/2 + 1) -n 


L ■ n 


jj add's 


L/2-n 


(L- 1) -n 



Table 1 . The operation cost of filtering a set of n points 
using the pair of filters if^, by 2 different techniques. 



^The additition of 1 to L/2 results from normalization by scalar- multiplication 
by Rl/2- 
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We note that by Theorem ^]T| the inverse wavelet transform has 



a decomposition given by reversing the order of the operators W2{aj) 
and replacing each by its inverse. 

3. Factorization of the biorthogonal wavelet transform. 

Let H'l, G| and be pairs of biorthogonal FIR analysis filters 

and synthesis filters, respectively. We make no restrictions on them, 
except that they satisfy the relations given in ( |1.13| ). 



We want to find a factoriztion of G„(if^, G^^) like we did in the 
case where the filters were orthogonal. This may seem slightly more 
complicated because of two major differences: 

• The filter lengths L, L of the lowpass and highpass filter need 
not be even numbers, moreover they will in general not be equal: 
L — L may be any number. 

• The highpass filter is not the alternating fiip of the lowpass 
filter H^, it is the alternating fiip of the synthesis lowpass filter 

But as we will see below, these differences are not crucial to the 
main idea in the method developed in the orthogonal case, and the 
relations ( 1.12|) , (|1.13|) will provide us with all that we need to make a 



modified version of the method work in this case. We simply have to 
remember that when dealing with the biorthogonal wavelet transform, 
we have to work with two "Multiresolution Analysises" (MRA's) which 
are duals of each other: One for analysis, the other for synthesis. Thus, 
to each linear operation/mapping we choose to perform on the analysis 
MRA, there corresponds a linear dual mapping on the synthesis MRA 
and vice versa. Having this in mind, we proceed to solve our problem. 

Now that we have the procedure from the orthonormal case to guide 
us, we begin with the definition of a linear map that reduces the length 
of the analysis lowpass filter H^. Since this will uniquely define some 
dual map on the synthesis lowpass filter, we see that we could equally 
well define our mapping on . It will be more practical to start with 

the shortest filter. We assume L > L. Since the filter if^ is not double- 
shift orthonormal to itself, F2{a) will only reduce this filter in one end 
for some a. Therefore, we may as well replace F2{a) by -^2(0^) to save 
useless operations. We define the first step analogous to 



Hi ^ F2{a,)Hl (2.15) 



where we claim 



Hl{0) ^ 0, which implies ai = -^y^. (2.16) 
From the "double shift biorthogonality relation" ( |1.12| ), we get 
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6o,k = {HlS{2k)Hl 

'F^\ai)F2{ai)HlS{2k)H 
'h{ai)HlF,-\a,)S{2k)Hiy 
Using relation ( |1.12| ) once again, we observe that 



= HliO)HliL-2) + HlmlCL-l] 
HliL-l) _ HliO) 

=^ F,-\a,):Hl{L- 1)^0. 
Figure illustrates this fact for the case L = 3, L = 5. 



(2.17) 



(2.18) 






• 


m m 


• 


• 


• 


• • 


• 

^1(4) 





• 


• • 

A(ai) 








9 


• 






Figure 3. Illustration of F2{ai)Hl and F2 \ai)H^. 
Then, using relation ( [L.13| ) and the duality of F2{ai) to F2~*(a;i) we 



get 



F2(ai)G|(0) = i-iy Fi\a,)Hl{L - 1) = 0. (2.19) 

Now we have to be careful in defining the next steps. We look 
for a sequence of linear non-singular operators of type F2(')? -^2(')) ^^e 
composition of which reduces the pair iJ^, G\ to a pair of 1-point fil- 
ters. By the relations (p.l7|) , ( p.l8|) , (|2.19|) , we see that this goal is 
reached once we have a sequence of pairs of non-singular dual linear 
maps that reduces the pair of filters H'l, if| to a pair of biorthogonal 
1-point lowpass filters. We will show that such a sequence exists by 
construction. 
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Using the operators -^2(0;), each step in this construction will reduce 
the filter in one end only. We have to decide which end to reduce 
at step k such that at every step j , I < k < j < max{L, L), 



supp{Hl^^) n supp{Hl(^^)) ^ 0, (2.20) 

and such that after the last step we are left with two dual 1-point filters. 

Since dual filters always have common center, we see that condi- 
tion ( |2.20| ) forces a symmetric reduction of the filter lengths around the 
common center. As is easily checked out, this is achieved by alternat- 
ingly making use of suitable operators Frn{c(k) , F^i^k+i) , m>2, 1 < 
k < max{L, L) — 1, together with suitable shifts. By the definition of 



the highpass filter given in (|1.13| ), we see that the centers of the 



lowpass and highpass filters do not share support, the highpass center 
is delayed 1 time step with regard to the lowpass center. Thus, the 
symmetric reduction of the highpass filter is delayed 1 time step with 
regard to the symmetric reduction of the lowpass filter. 

We note one possible problem: If at some step j, I < j < max{L, L) 
the filter H'^_^ has internal zero coefficients, we will not be able to map 
the outer (with regard to the center of the filter) neighbour coefficient 
to zero using any F2{c() or -^2(0^)- But it is easy to see that by replacing 
F2{a) by -^3(0;) or -F|(a) by F^{a) , as defined in (|2.3|), the desired zero- 



mapping of the neighbour coefficient can be achieved in both H^_^ and 
H^jy In general, if there are m consecutive internal zero coefficients, 

one should use Fm+2(«), -F4+2('^)- 

There is really not much left to prove, except to clear up some 
details. We separate the biorthogonal filters into 3 classes, and give 
the form of our factorization of the wavelet transform in each of these 
classes. 

3.1. The general odd-odd case. Here L 7^ L are both odd num- 
bers. To simplify notation we define 



^_ L + L 
2 



yymik){ak) = S r^t \Q{ k 



^„^{fc)(afc)5(-[|J) if A;, 1 < A; < J, is odd. 
^Lik)MSi-l) if A;, 1 < A; < J, is even. 

(2.21) 



Then we can state the following result: 

Theorem 2.2. Given a pair of biorthogonal FIR filters H^, G| of 
odd lengths L,L, there exists a sequence of integers {fTi{j)}j^-^ C N — 
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{1} such that the operator <dn{Hi^,GY) defined in ^2. 1\ ) may he decom- 
posedasQn{HlG\) = S{-[J/2\)W-lj^{aj)W^^^^^ 
with Wm(j){,(y-j) , ^ ^ j ^ J (ind J defined as in ( \2.21\ ). 



Proof. The theorem was proved by construction. The argumen- 
tation below is just a formahzation of this construction. 

cl = {^*2crHl){k) = J2^{j)H{j-2k) = {S{2k)Hl^) 

= ■ ■■W-l^{aj)W^^j^{aj) ■ ■■W^(,){a,)S{2k)Hl^) 

= {WmiJ){aj) ■ ■ ■ Wmii){a,)S{2k)Hl W-^^iaj) ■ ■ ■ #-(\)(ai)x) 

= {6{j, 2k - J/2), W-l^{aj) ■ ■ ■ W^~(\)(ai)x) 

= ■ --W-l^i^iM^k - LJ/2J). (2.22) 

4 = {^*2aGl){k)=Y,^ij)GlU-2k) = {Si2k)Gl^) 

■ --W-l^iajW^^j^iaj) ■ ■ ■ W^^(i)(«i)5(2A;)G^,x 
W^iJ)iaj) ■ ■ (ai)S(2fc)G|,#-f^)K) ■ ■ ■^^^(aOx 

[6{j, 2k + 1- J/2), W-l^iaj) ■ ■ ■ #-(\)(ai)x) 
= W-l^{aj)---W~l^{aM2k+l-[J/2\). (2.23) 

It remains to prove that the orthogonal length-reduction of the filters 
H'[, iJ| takes steps. But this is easy to see: We have 



by relation (|1.12| ). If L — L = 2(2A;+ 1), A; > 0, then by considering rela- 
tion (|1.12|) we deduce that 2k of the L — L + 1 coefficients in iJ|_^^j^ are 



0, and that these 0-coefficients are distributed symmetrically around 
the center-coefficient. Thus the total number of steps will be 

L-L L+L 

L-l + 2k + 2 = L-l + + 1 = . 

2 2 

By a similar argument, if L — L = Ak, k > 1, then the total number of 
steps will be ^ - 1. 

Thus, the number J as defined in ( |2.21| ) is at least an upper bound 
for the total number of steps. 

□ 
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A schematic picture of the algorithm that results from the factor- 
ization of 0„(if|,G'|) in this case of odd length filters is shown in 
Figure ^ below. At level k in the figure, each slanting line A connected 
to a vertical line B means multiplication of the number in the point 
of origin of A by the factor a^, followed by addition to the number in 
the point of origin of B, while single vertical lines just mean the iden- 
tity. This is easily seen to be the pairwise mapping induced by either 



or F2 (ofc) 



r/i A 

Cq Uq c-^ 



x(0) x(l) 



"n/2-1 



u 



x(n- 1) 



Figure 4. Filtering-diagram for a conjugate pair of 
biorthogonal filters of length 5 and 3. 



The "tilde" c's and d's at the edges of the filtering-diagram above 
can be obtained by mirroring x about its endpoints or periodizing. 
Counting operations, we see that we have to do 1 multiplication and 1 
addition pr. pair of points for each Wm{j){(^j)- Writing 



1 

n„(iJiG|)^n^4)("^)' (2.24) 

i=J 

then by construction n„ = >S'([^J)6„, but it has a more efficient im- 
plementation as shown in Table ^ 



Operation 


n„ 




jj mult's 


4 




% add's 


L+L 
4 





Table 2. The operation cost of filtering a set of n points 
using a pair of biorthogonal odd-odd filters H'l, by 2 
different techniques. 



^The cost result shown in Tabic ^ is only an upper bound if the H'l, G^^ have 
internal zero coefficients. 
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We note that by Theorem the inverse wavelet transform has a 



decomposition given by reversing the order of the operators W2 (ctj) 
and replacing each by its inverse. 

3.2. The symmetric odd-odd case. There is only one crucial 
observation to be made here, then the results follow immediately from 
our work in the previous section. This is best illustrated by an example. 
Figure ^ shows the first 2 steps in reducing a symmetric filter of length 
7 to a 1-point filter. 






Figure 5. 2-step reduction of a symmetric filter of odd 
length using some F2(ai), -p2('^2)'S'(— 1). 

By the definition of our length reducing maps in the section above, 
we have a2k-i = C(2k, 1 < ^ < max{L, L)/2, because the maps 
F2{c(2k-i)S{—{k — l)), F2{a2k)S{—k) preserve symmetry in pairs. Thus 
we see that replacing the operator F2{a2k)S{—k)F2{a2k~i)S{—{k — 
< k < max{L,L)/2 by the operator F3{ak)S{-{k - 1)),1 < 
k < max{L^L)/2 as defined in (|2.4| ), will improve efficiency: F^{ak) '■ 

Tjb I Tjb 

^L-2k ^ ^L-2k-2- 

Using similar notation to the previous section, we define 



W^{ak) = F,{ak)S{-{k - 1)). (2.25) 
Then we have the following refinement of Theorem 



Corollary 2.1. Given a pair of biorthogonal symmetric FIR fil- 



ters H^,G'^^ of odd lengths L,L, Qn{FI^,G^^) defined in ( ^.^ may be 
decomposed as 

Qn{HlG\) = S{-\'^-\)W^\aj)W.^\aj.{) ■ ■■W.^\ai), with W^,), 
1 < j < J and J defined as in ^2.23i ). 

Proof. The only thing left to prove is that the total number of 
steps equals J as defined in ( |2.25| ). But this follows easily by the same 



argumentation as in the proof of Theorem |2l2. 



□ 

A schematic picture of the algorithm that results from the factor- 
ization of Bn(-ffi,G|-) in this case of symmetric odd length filters is 
shown in Figure |^ below with the same use of symbols as in Figure ^ 
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x(0) x(n - 1) 



Figure 6. Filtering diagram for a symmetric pair of 
conjugate filters of length 7 and 5. 

The "tilde" c's and d's at the edges of the filtering-diagram above 
can be obtained by mirroring x about its endpoints or periodizing. 
Counting operations, we see that we have to do 1 multiplication and 2 
additions pr. pair of points for each Wm{j){aj)- We write 

1 

n„(i7^,G|) = n^i)(«^)- (2.26) 

By construction n„ = >S'([^] — 1)6„. Iln has a more efficient im- 
plementation as shown in Table ^. We note that if the filters H, G 
have inner zero-coefficients, we simply replace -^3(0) by the composi- 
tion some m > 3, and it is easy to see that this does not 
affect the total number of additions and multiplications. f\ 



Operation 






tl mult's 






jj add's 







Table 3. The cost of filtering a set of n points using 
a pair of biorthogonal symmetric odd-odd filters if^, G*! 
by 2 different techniques. 



3.3. The even-odd case. We assume without loss of generality 
that L is odd and L is even. We observe that there is only one n such 
that Hl{n)H\{n) ^ 0. 

Now, define 



■^The cost result shown in Table || is only an upper bound if i7^,G^ have 
internal zero coefficients. 
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^_ L + L-l 
2 

W'„(j)(aj) as in (2.21) (2.27) 

(2.28) 

Then we have the result: 

Corollary 2.2. Given a pair of biorthogonal FIR filters H\-,G^^ 
of odd- even lengths L, L, there exists a sequence of integers {Tn{j)}'j^^ C 
N — {1} such that the operator 0„,(if|^, G|) defined in ( ]^. i| j may be de- 
composed as 

BniHlGl) = 5(1- LiJ)l^2"*(«j)W^-(Vi)(«^-i)---W^™fi)("i)' ^^^^ 
W^m{j)(«i); ^ ^ j ^ J (ind J defined as in 



Proof. We only need to show that the total number of steps equals 
J as defined in ( p.28|) . We have 

Now, setting L — L = 2A; + 1, A; > 0, we observe that Hl_^^^ has only 



+ 1 nonzero coefficients because of relation (|1.12|) . This yields 



L-L+1 L+L-l 
J = L-l + k + l = L-1 + = . 



□ 



Counting operations, we see that we have to do 1 multiplication 
and 1 addition pr. pair of points for each Wm{j){aj), 1 < j < J- 

1 

n„(ifi,G|)^n^4)("^)' (2.29) 

j=J 

then by construction ^(l — [^J)n„ = Qn- n„ has a more efficient 
implementation as shown in Table ^. 

3.4. The even-even case. Here, both L, L are even numbers. 
There is only one important difference in this case: The last step will 
not reduce the dual filter from a 2-point filter to a 1-point filter, as 
illustrated in Figure 0. 

This fiaw is easily repaired by replacing ^2(0^) by F2{ak,P) in the 
last step k, where P is chosen such that F^'^^ak, yields a 1-point dual 
filter. Using similar notation to the previous sections, we define 

"^The cost result shown in Table ^ is only an upper bound if i7^,G^ have 
internal zero coefficients. 
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Operation 






jj mult's 


L+L-l 
4 " 


¥^ 


tl add's 


4 





Table 4. The cost of filtering a set of n points using a 
pair of biorthogonal even-odd filters G| by 2 different 
techniques. 







1/S 



Figure 7. The result of reducing a pair of dual 
biorthogonal filters of length 4 to 1 point. 



j_ L + L 
^" 2 ' 

r F2(«j,/?)^(-(J-l)/2) if A; = J. 
= <^ FUak)S{- LfJ ) if A;, 1 < A; < J - 1, is odd. 

[ F^{ak)S{-[!^\) ifk,l<k<J,is even. 

(2.30) 

Then we have the following result: 

Corollary 2.3. Given a pair of biorthogonal FIR filters 
of even lengths L,L, there exists a sequence of integers {fn{j)}j^i C 
N — {1} such that the operator Qn{H^,G^i) defined in (\2.1\ ) may he 
decomposed as 

Qn{HiG\) = 5(1- LiJ)i^2"*(«j)w^-(Vi)(«^-i)---w^™;i)(«i). y^^ih 

Wm{j){oij), I < j < J and J defined as in (^-SQ ). 

Counting operations, we see that we have to do 1 multiplication 
and 1 addition pr. pair of points for each Wm{j){aj), I < j < J, and 2 
multiplications and 2 additions pr. pair of points when implementing 
Wj{aj). Writing 

1 

n„(i7iG|) = n^4)(«^)' (2.31) 
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then by construction ^(l — [■fj)n„ = 9„. n„ has a more efficient 
implementation as shown in Table |. f\ 



Operation 






jj mult's 


L+L+6 
4 


¥^ 


tl add's 


L+L+2 
4 





Table 5. The cost of ffitering a set of n points using a 
pair of biorthogonal even-even filters H'l, G| by 2 differ- 
ent techniques. 

Remark. Since orthonormal FIR ffiters are a special case of the filters 
discussed in this section, we see that it is possible to reformulate the 
theory for orthonormal filters by exclusively using elements in SL2(R). 
Furthermore, it is possible to achieve the same operation cost result as 
before. 

We note that in is shown fast algorithms for computing wavelet 
coefficients with similar operation-counts to ours using a different ap- 
proach. 



^The cost result shown in Table || is only an upper bound if i7^,G^ have 
internal zero coefficients. 
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CHAPTER 3 



Extension of Results From Dimension 1 to Higher 

Dimensions. 



When considering tensor wavelet bases in dimension 2 or 3, one 
will face the problem of filtering a set of n x n or n x n x n points, 
respectively, in an effective way, where n is some power of 2. Assuming 
identical MRA's in all dimensions and using notation as in the previous 
chapter, we look for efficient implementations of 6„(iJ, G) ® 0„(//, G), 
Qn{H, G) ® Qn{H, G) <S) Qn{H, G). Proceeding straightforward, we can 
use the factorization of the wavelet transform in dimension 1 separately 
in each of the 2 or 3 dimensions, and thus achieve the same operation 
cost result as before: Reducing the operation count roughly by half. 
However, in both the orthonormal and biorthogonal case it is possible 
to obtain substantially better results. 



2.1. The orthonormal case. Since our factored wavelet trans- 
form in dimension 1 only works with 2 points at a time, and tensor 
products commute, we can restrict to a square consisting of a set of 
2x2 points at a time, as illustrated in Figure ^ below. 



1. Introduction. 



2. Dimension 2. 



(2) 




I 2 I 



2 



2 




I 2 I 



2 



2 




xo,i xi,i 



2 



2 




Xo,0 Xi 



(1) 



Figure 1. Filtering-diagram for the discrete wavelet 
transform in dimension 2. 
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The Figure D illustrates the division of the point-set into point-squares. 
Each point-square represents the map induced by restricting some 
F2{a) <S) F2{a) to this point set, more specificly each line connecting 
a pair of points represents the map induced by restricting some F2{a) 
to this pair of points. The squares will alternatingly be in the posi- 
tion marked 1 and 2 as A; in F2{ak) ® F2{ak) varies. It will turn out 
favourably to use the matrix M2{a) on its trigonometric form, that is 



^ ' \^ sma cos a J ^ ' 

Using this form and restricting to a point-square, we write 

F2{a)®F2{a)\,^2 ■ { ] - ( j • (3-2) 

i xo,i xi,i j [ yo,i yi,i j ^ 

Computing on a square yields 



f yo,o \ ^( cos^a(xo,o -xi,i) - sinacosa(xo,i + xi,o) +xi,i 

V yo,i / V ~" ^^^^ "(^0,1 + xi,o) + sin a cos a(xo,o - xi^i) + xo,i 

( ^I'O ^ = ( -sin2a(xo,i + xi,o) +sinacosa(xo,o -xi^i) +xi,o 

V yi,i J \ - cos^ a(xo,o - xi,i) + sin a cos a(xo,i + xi,o) + xo,o 



Writing 



we get 



d = cos^q;(xo,o — Xi_i) — sina;cosQ;(xo,i + Xi^o), 



yo,o = d + xi^i 

yo,i = cot a ■ d + xo,i 

yi,o = cot a ■ d + xi^o 

yi,i = -d + xo,o. 

Combining these relations, we see that ^2(0;) ® ^2(0) [2x2 can be im- 
plemented by 3 multiplications and 7 additions. Using notation as in 
dimension 1, we write 



1 

^niHlGi)^ n (W^2K)®W^2K))l2x2> (3-3) 

j=L/2 
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Then by construction H„ = S{L/2 — 1)G„ S{L/2 — l)9n, but S„ has 
a more efficient implementation, as shown in Table |l|. 



Operation 


' — "n 


n„ (g) n„ 


On ® 0„ 


jj mult's 




(L + l)n2 




jj add's 






(2L - 2)n^ 



Table 1 . The operation-cost of filtering a set of n x n 
points using the pair of ffiters Hf, Cl by 3 different tech- 
niques. 



2.2. The biorthogonal case. Arguing as in the last section and 
using the same notation as in dimension 1, we write 



^2 







_ / yo,o 


yi,o 1 


2x2 


\ ^0,1 Xi 1 J 


I xo,i 


yi,i J 



(3.4) 



a 



3x3 




^0,0 yi,0 ^2,0 

yo,i yi,i y2,i 

Xo,2 yi,2 X2,2 



(3.5) 



We see that for even-even filters and general odd-odd filters we can 
restrict to a set of 2 x 2 points at a time, while for symmetric odd-odd 
ffiters we restrict to a set of 3 x 3 points at a time. This "factorization" 
of operations allows us to implement the biorthogonal wavelet trans- 
form more efficiently in dimension 2, as illustrated in Figure 0. Each 
arrow-line in the figure means multiplication of the point of origin of 
the arrow by the number a, followed by addition of the product to the 
point of termination of the arrow line. 



To implement F2 



a) efficiently, we compute on a square 



of 2 X 2 points as shown in Figure |^. Writing 



yo,o = xo,o + axo,i 

yo,i = xo,i 

yi,i = xi^i + axo,i 

yi,i = xi,o + a(yi,i +xo,o), 



-•^The addition of 1 to L in the multiply-cost for n„ (g) n„ results from scalar 



multiplication by R 



L/2- 
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(2) 



Xo,2 
• 



Xo,l, 



Xl,2 
« 



X2,2 



Xl,l 



X2,l 




Xll 



Xo,0 



Xl,0 



X2,0 



Xo,0 



Xl,0 



(1) 



Figure 2. The point-squares for symmetric and non- 
symmetric biorthogonal filters, respectively. 



we see that we can implement F2 {a) ^ ± 2 



a, 



2x2 



by 2 multiplica- 



tions and 4 additions. 

Similarly, to implement F^^{a) (8> F^^{a) efficiently, we write 



yo,o 


= Xq^o 


yo,i 


= Xo,l 


yo,2 


= Xo,2 


yi,2 


= Xi^2 


y2,2 


= X2,2 


y2,o 


= X2,o 


yi,o 


= Xi^o + "Xo,o + «X2,o 


y2,i 


= X2,l + aX2,2 + ax2,o 


yi,i 


= xi^i + a(yi,o + ax2,2 + X2,i 



Thus, we see that F. (g) Fo *l 



a. 



3x3 



can be implemented by 2 mul- 



tiplications and 8 additions, assuming that we have stored the numbers 
ax2,2, axo,o, 0x2,2 + xi,2 from previous steps in the implementation of 
the algorithm. 

Using notation as in ( p.21| ), in the general odd-odd case we define 



J](W^2K-)®W^2(«,)) 

and thus get S„ = S([J/2j)e„ ® 5([J/2j)e„. 



2x2 



(3.6) 
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In the odd-even case we define 



2x2 



where we use the notation defined in ( |2.28| ). We get 
In the even-even case we define 



(3.7) 



2x2 



j=J 

where we use the notation defined in (|2.3CI|) . We get 

i)e„®5(L|j -i)e„. 

Finally, in the symmetric odd-odd case we define 



(3. 



3x3 



(3.9) 



where we use the notation in ( p.25|) . This yields S„ = S'([-^])6„ ® 

Thus, we get the operation-cost results shown in Table ||, Table H, 
Table |. g 



Operation 


' — 'n 




e„ (g) e„ 


tl mult's 


L+L 2 
4 




(L + L)n^ 


tl add's 


^+^n2 




{L + L- 2)n2 



Table 2. The operation cost of filtering a set of n x n 
points using a pair of biorthogonal odd-odd filters H'l, 
by 3 different techniques. 



3. Dimension 3. 

It is straightforward to see that we can group the operations of 
F2{a) (g) F2{a) ® ^2(0), F2{a) O ^2(0) O ^2(0) into cubes of 2 x 2 x 2 
points, and F3 (a) 0^3(0) (8)^3(0;) into cubes of 3 x 3 x 3 points. This is 
illustrated in Figure ^. To make things easy, we will use some "loose" 
geometrical terminology. 

^The operation-counts shown in Tabic |[ Tabic ||, Table U, Table are only 
upper bounds if H^, G| have inner zero coefficients. 
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Operation 


' — "n 


n„ (g) n„ 


e„ ® ©n 


ji mult's 


L+L-l 2 
4 " 




(L + L)n^ 


jj add's 






{L + L- 2)^2 



Table 3. The operation cost of filtering a set of n x 
n points using a pair of biorthogonal odd-even filters 
H'I, G*! by 3 different techniques. 



Operation 




n„(g)n„ 


©„®©n 


tl mult's 


L+L+10 2 
4 


L+L+4^2 


(L + L)ri2 


jj add's 


L+L+2^2 


L+L+2^2 


{L + L- 2)n^ 



Table 4. The operation cost of filtering a set of n x 
n points using a pair of biorthogonal even-even filters 
H^,G'f by 3 different techniques. 



Operation 






©n ® ©n 


tl mult's 




L^Jn2 


(L + 


jl add's 




2L^Jn2 


(L + L - 2)^2 



Table 5. The operation cost of filtering a set of n x n 
points using a pair of biorthogonal symmetric odd-odd 
filters H^, G| by 3 different techniques. 



3.1. The orthonormal case. We did not succeed in finding any 
symmetries in the trigonometric expressions in the 3-dimensional case. 
The cheapest arrangement of operations will therefore be to apply 
W2{a) W2{a) |2x2 to 2 "opposing planes" of 2 x 2 points in the 2x2x2 
cube, and then W2{<y) to the 4 "lines" that connect these 2 planes. 
Thus, one has to do 14 multiplications and 22 additions pr. cube, not 
counting normalization. We denote this "factored" wavelet transform 
in 3 dimensions on nxnxn points by fi„(Lr^, G^)- Like before, we have 
Qn = S{L/2 - 1)©„ ® S{L/2 - l)Qn ® S{L/2 - 1)©^„. The operation 
cost-result is shown in Table ^. 

3.2. The biorthogonal case. By rearranging the arithmetic op- 
erations, we find we have to do 6 multiplications and 12 additions pr. 
cube in the even-even/odd-even and non-symmetric odd-odd case. In 

•^The addition of 1 to |L and |L is because of normalization by R\i2- 
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Figure 3. The cube for biorthogonal odd-odd/even- 
even filters on the left and for symmetric odd-odd filters 
on the right. 



Operation 


^71 






jl mult's 








tl add's 






3(L - l)n^ 



Table 6. The cost of filtering a set of x n x n points 
using the pair of filters G^, by 3 different techniques. 



the symmetric odd-odd case we found an operation arrangement yield- 
ing 6 multiplications and 24 additions pr. cube. All these operation 
counts are easily verified by considering Figure 0. Denoting in each 
case the factored biorthogonal discrete wavelet transform on n x n x n 
points by we get that up to shifts, f2„ equals 9„®9„®0n- 

The cost-results are shown in Table ^ Table ^, Table ^ and Table |10. 



Operation 








tl mult's 


1{L + L)n^ 


\{L + L)n^ 


|(L + L)n3 


tl add's 


1{L + L)n^ 


^{L + L)n^ 


|(L + L-2)n3 



Table 7. The cost of filtering a set of x n x points 
using a pair of odd-odd filters H^, G^, by 3 different tech- 
niques. 



^If iJ^,G| has inner zero-coefficients, the operation counts in Table Table 
^, Table ^ and Tabic |l^ arc only upper bounds. 
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Operation 




n„ n„ ® n„ 


©n ® On ® 0„ 


mult's 




f(L + L-iy 


|(L + L)n3 


add's 


f(L + L-l)n3 


f(L + L-l)n3 


|(L + L-2)n3 


Table 8. The cost of filtering a set of n x n x n points 
using a pair of odd-even filters H^, G^^, by 3 different 
techniques. 


Operation 




n„ (g) n„ (g) n„ 


0„ (g 0„ ® 0„ 


jj mult's 


i(L + L + M)n3 


f(L + L + ^)n3 


f(L + L)n3 


tt add's 


1{L + LW 


f(L + L + 2)n3 


f(L + L-2)n3 



Table 9. The cost of filtering a set of n x n x n points 
using a pair of even-even filters by 3 different 

techniques. 



Operation 






©n (g ©n ® ©n 


mult's 








add's 




3L^^Jn=^ 


|(L + L-2)n=^ 



Table 10. The cost of filtering a set of n x n x n points 
using a the pair of symmetric odd-odd filters i?^, by 
3 different techniques. 



CHAPTER 4 



The 5'02(R)-Method in Computing Filters and 

Edge Maps 

In this chapter we exploit our 2-point factorization of the orthonor- 
mal discrete wavelet transform in the last chapter to compute filters by 
solving non-linear equations in several variables numerically. We also 
try to improve on the preservation of regularity in the discrete wavelet 
transform by constructing some (non-orthogonal) edge maps. 

1. Computation of orthonormal filters 

We start with the construction of an orthonormal filter of length 
4 with a maximum number of vanishing moments compatible with its 
support width. That is, we want ai, 0:2 so that 



F2(«2)5(l)F2(ai)l ^ 6^,k 
F2(a2)5'(l)F2(ai)n ^ 6i^k 

This is simple enough to do by handcalculation. 

^ aix(O) + x(l) — a;2x(2) + aiQ;2x(3) 

1 _ aia;2x(0) + a2x(l) + x(2) - a;ix(3) 



(4.1) 



(4.2) 



(4.3) 



We insert the vanishing moment conditions from (|4.1| ) and get the 
equation-set 



aia2 + 02 ~ «i + 1 = 
a2 — 3ai + 2 = 



Solving yields 



G {(I/V3, - 2), (-I/V3, 2 + V3)} 

We find the lowpass filter by inserting the solutions in (|4.2| ). 

The two solution-sets are mirror-images of each other, it is the filter 
shown in Table |1| that corresponds to our choice of rotation-matrix 

M2(«fc). 
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n 




Ctn 





0.482962913 




1 


0.836516304 




2 


0.224143868 


\/3- 2 


3 


-0.129409523 





Table 1. The filter H4 and the entries in the corre- 
sponding M2{an)- 



To find orthonormal filters of length L, the above method leads 
to the wearying task of solving ^ nonlinear equations in ^ variables. 
However, by use of some Maple-routines we succeeded in generating 
Hq and some orthogonal filter Hg corresponding to a wavelet with 4 
vanishing moments. 



n 


Hg(n) 







.03222310057 




1 


- .0126039675 


- 2.556583915 


2 


- .09921954341 


- .1434214911 


3 


.2978577953 


.7958755204 


4 


.8037387510 


- 2.351285662 


5 


.4976186668 




6 


- .02963552763 




7 


- .07576571472 





Table 2. The computed filter Hg and the entries in the 
corresponding M2(afc). 



Remark: The filter in Table is not the one given in |T| , but these 
filters are not uniquely determined by claiming a maximum number of 
vanishing moments only. 

Now we turn to coifiets. Claiming vanishing moments on the scaling 
function leads to the following observation: 

If H{n) ■ = for V p, 1 < p < M < 1, 

n=0 ^ 

where M is odd, then 

L-l 

^/7(n) ■n*'+^ =0. (4.4) 

n=0 

This means that it suffices to impose odd vanishing moments on 
0. Counting degrees of freedom in the filtering-diagram we find that 
for ^2°*'^'^*, if the highpass-filter has M vanishing moments, then the 
lowpassfilter has 2(-| — M) vanishing moments. The proof of this is 
given below. Maple were able to solve the nonlinear equations to give 
jjcotfiet some ^ where we used 3 out 4 degrees of freedom 



1. COMPUTATION OF ORTHONORMAL FILTERS 
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to provide the highpass filter with 3 vanishing moments, yielding 2 
vanishing moments on the scaling function. 



n 






-3 


- .07342587213 




-2 


- .03314563036 




-1 


.4854426593 







.8065436723 




1 


.3100524696 


2.215250436 


2 


- .09943689110 


.1504720765 


3 


- .01496247550 


- .7384168123 


4 


.03314563037 


- .451416230 



Table 3 . The computed filter Hg°^^ ^ and the entries 
in the corresponding M^ipin)- The highpass filter has 3 
vanishing moments, the lowpass filter has 2 vanishing 
moments. 



Proof of ( Ml) 
We have 



= y x^(f){x)dx ,Vp, l<p<M<^-l 



5=0 



0, for V j9, 1 < p < M < - - 1. 



Fourier transforming equation (|1.1| ) we get 



0(0 = mo(e/2)0(e/2) 

where mo is the trigonometric polynomial that is the Fourier transform 
of the filter i?™*'^- Differentiation yields 



= 0'(O) = im[,(O)0(O) + ^mo(O)0'(O) = ^ 
Repeating this argument M times we get 



d'P 
d^ 



rriQ 



0, for V p, 1 < p < M, 



5=0 



thus we conclude that itlq must be on the form 



mo(0 = ! + (!- e-^«)''^'P(0 



(4.5) 
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for some trigonometric polynomial P. Assuming that the coefficients 
i?™*'^(n) of mo are real, we have 



moi-O = ^ moiO -i-^ mo(0 = moiO, 

which implies 

Expanding mo in a Taylorseries around the origin 

mo(0 = 1 + ai^ + a2^^ + ■ ■ ■ + a„C + ■ ■ ■ ,where = ° / 
we get 

imo(oi' = (i + aie + a2e' + 03e^ + ---)(i-«i^ + «2e'-03e^ + --- 

— J- + CaZ+i^ + + ■ ■ ■ + Cm+tiC, + " " " 

for some sequence {cM+i}iZi because of ( |4.5| ). By computing the above 
product of the two series, we get the relations 

Co = 1 , C2fc+l = 0, 

fc-i 

C2k = i-'^)^al + 2 ^(-l)*a2fc-iai, where Uq = 1. (4.6) 

i=0 

Moreover 



\akk\\ = |m;,'^(0) 



(4.7) 

Now, claiming the Ck = for 1 < < M and using ( [4.6|) , ([4.7|) , the 
result follows. Indeed, one has 

^-1 T 

H'^'^{n)nP = for V p, 1 < p < M < - - 1, 

n=0 

where M is odd, 

L-l L-l 
n=0 n=0 

□ 



2. WAYS OF DEALING WITH THE EDGE PROBLEM 
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2. Ways of dealing with the edge problem 



We look for ways of "smoothing" the edges of a vector in R". That 
is, we look for operators acting on the edge coefficients that to some 
degree possess the following properties: 

• PI: Linearity and orthogonality. 

• P2: In some strict sense, "nice functions" should map to "nice 
functions" . 

• P3: The map should ensure some vanishing moments on the c/'s. 

We have to make clear the meaning of P2: Considering the coeffi- 
cients Cfc at some level j, we see from 



that a polynomial maps to a polynomial of the same degree under 
QniHiiGb)- Depending on the number of vanishing moments on the 
corresponding if), the rf^'s will be zero. 

2.1. The method of periodizing. This method simply peri- 
odizes the vector around its edges. This yields orthogonality of the 
transform, but only preserves regularity up to constant functions, and 
we get only 1 vanishing moment on the d^s. 

2.2. The method of mirroring. One way to overcome the edge 
problem is to mirror the input vector about each of its endpoints and 
then apply the filter to the mirrored points. 

The method of mirroring does not lead to an orthogonal transform, 
but it preserves regularity up to continuity, and yields 1 vanishing 
moment on the d^s. 

2.3. The method of edge matrices. Another approach is to try 
to construct a matrix map acting directly on the edgepoints such that 
the properties PI, P2 and PS are satisfied to some extent. We will 
do this construction by claiming some vanishing moments on the d^s 
and claiming some degree of preservation of regularity of the trans- 
form at edges by claiming monomials up to some degree M map- 
ping continuously to polynomials of degree M. That is, we com- 
pute the polynomials Pj given in ( [4.8[ ) for < j < M, and impose 

Pj ,0 < j < k{L) < M at edges, L is the length of lowpass filter. 




Puik) , Pm a polynomial of degree M, 



(4.8) 
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For a orthonormal FIR filter Hl we will construct two — 1) x — 1) 
matrices Eh^a, EHL,r, each acting on the ordered set e^, of ^ — 1 
edgepoints on the lefthand and righthandside, respectively, to give the 
"tilde" coefficients. For example 

E//io,« ■ 6/ (co, do, ci, di) 

E//i(,,r : ^ (Cn/2-2) dn/2-2, C„/2-l, dn/2-l) 

This is illustrated in Figure |1|. 



Co (io Ci di c\ 




Xo Xi X2 X3 



Figure 1. 

As for fulfilling the above listed properties, this edge map is of course 
linear, but not orthogonal. We have sacrificed orthogonality to be able 
to achieve P2 and P3. Some edge-matrices are shown in the appendix. 

2.4. Comparing the methods. The profits of the edge matrices 
are obvious: 

• Good: Preservation of regularity to a higher degree than by 
simply mirroring at edges. 

The drawbacks are equally obvious: 

• Bad: This edgemap is very far from orthogonal in the sense that 
the operator norm of the map turns out to be large, and it seems 
to grow with increasing degree of preservation of regularity as 
can be seen in Table ^ and Table ^ 

To be able to compare the methods of mirroring and edge matrices 
in this non-orthogonality respect, we proceed as follows: 

Given the vector x G R" and some orthogonal FIR filter Hi, we 
compute the c^'s and the c/^'s, getting 



K{k) 

Ck = Y^ 7(i)x(i) 

j=0 



(4.9) 
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when we use edge matrices, and 
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j(fe) 

i=0 



(4.10) 



when we mirror x. The expressions for the d^s are similar. We use 
Cauchy-Schwartz on (f4.9|) and ( [4.10| ). On the left edge we set 



Qiick) 



PliCk) 



K{k) 
i=0 



j(fc) 

E 

i=0 



1/2 



1/2 



On the right edge, the bounds Qr, Pr are defined similarly. Now, 
we may compare the two methods in the non-orthogonality respect. 
Results corresponding to the edge matrices shown in the appendix are 
given in Table ^ and Table |^. 



P bounds on 
the \ck\, \dk\- 


L in Ht- 


6 


8 


10 


12 




.852 


.935 


1.05 


1.04 


Qiico) 


9.50 


21.8 


310 


757 




.989 


.907 


.973 


.963 


Qi{do) 


1.39 


2.14 


11.0 


2.02 


Pi{ci) 






1.04 


1.12 


Qi{ci) 






48.0 


103 


Pi(di) 




.994 


.999 


.996 


Qi{di) 




2.59 


1.58 


7.98 


P{d2) 








1.00 


QM2) 








1.10 


Pr(co) 


.852 


1.01 


1.05 


1.01 


(3r(co) 


1.01 


1.03 


1.01 


1.03 


Pr{do) 


.989 


.994 


.973 


.963 


Qr(rfo) 


.512 


.402 


.973 


3.13 


Pr(ci) 




.935 


1.04 


1.04 


(3r(ci) 




1.00 


1.00 


1.00 


Pr(rfl) 






.999 


.996 


Qr(dl) 






2.58 


5.03 


Pr(52) 








1.12 


Qr(c2) 








1.00 



Table 4. The bounds for the edge operators Ej^^d . 
using Daubechies shortest filters. 
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P bounds on 


L 


in m"'' . 


the Cfc|, |dfe|. 










6 


8 


12 


Pi{5o) 


1.13 


.949 


1.16 


Qi{cq) 


4.75 


2.43 


25.0 


Pi{do) 


.930 


1.15 


.909 


Qi{do) 


1.97 


10.4 


2.31 


Piici) 






.997 


Qi{ci) 






2.50 


Pid^) 




1.02 


.999 


Qi(di) 




1.54 


2.35 


P{d2) 






1.00 


QM2) 






2.31 


fr(co) 


1.13 


.907 


.997 


Qr(co) 


1.02 


1.97 


1.45 


Pr(do) 


.930 


1.02 


.999 


Qr(do) 


.705 


.908 


3.33 


Pr(£l) 




.949 


1.16 


Qr{cx) 




1.10 


1.03 


Pr(dl) 






1.00 


Qr{dl) 






19.5 


Pr{c2) 






.997 


Qr(c2) 






1.16 



Table 5. The P bounds for the different edge operators 
E^coi/ using Coiflet fihers. 

Table ^ and Table |^ clearly show the blowup of some of the edge- 
coefficients that results from using edge matrices. However, for the 
Daubechies filters the blowup effect almost exclusively affects the left- 
handside, while for the Coiflet filters the blowup effect does not seem 
to favour any side, and is on the whole more moderate. Anyway, this 
blowup will affect the numerical stability of the filtering-reconstruction 
procedure. 



CHAPTER 5 



Classification of Radar Signals Using Local 
Feature Extraction in the Space-Frequency Plane 

1. Formulation of the problem 

We consider the problem of separating two different distributions 
(classes) of electromagnetic radar signal sources from one another by 
doing a space-frequency analysis on the signals. Each distribution 
Dn{a) of signal sources will consist of a number of n coherent signal 
transmitters distributed randomly over n regularly spaced plane do- 
mains {fln,i{c()}i'=i with one and only one transmitter in each domain 
fl{n,i). The domains may overlap. 



f]„ .(a) = {(rcos^,rsin^) : 1 < r < 10, 2tt- <e < 27r- + a}. 

n n 

A picture of D^l-) and -D4(-) is shown in Figure |l], where the shaded 
areas are the {^3A}i=i and the {f^4,i}f=i. 



i"(inner) ^ 1 



r(outer) = 10 




r(inner) ^ 1 



r(outer) = 10 




D(4) 



D(3) 



Figure 1. The domains of D^^-) and D^^-). 

The analytical forms of the signals are solutions of the wave equa- 
tion in three dimensions: 



where 



- c^V^ip{^,t) = 0, 



(5.1) 



X2 



92x3 ■ 
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58 5. CLASSIFICATION OF RADAR SICNALS 

Making the ansatz ilj{'K,t) = 0(x)e~*'^*, we reduce the problem to 
the Helmholtz equation 

V^(/)(x) + fc20(x) = 0, k = -. (5.2) 

c 

The fundamental solution of this equation may be expressed in polar 
coordinates as 



0(r, 6) = (f)ir) = A , A is a constant. 

r 

Writing {pi}"^]^ for the locations of the signal sources in a signal 
source distribution Dn{-), we get the following analytical expression for 
the signal s„(x, t) corresponding to this distribution: 

s„(x,t) = e-^^^y -^^e''^\^-P< (5.3) 

We write 



X = {Rcos9, Rsm9), pi = (rj cos 6'j, sin 6^4). 
Now, we assume R is large compared to the r^'s and make an expansion: 



|x-pi| = R 



1 + 



rf — 2_Rr j (cos 9 cos 9i + sin 9 sin 9i 



i?2 



1/2 



R + — r i{cos 9 COS 9i + sin 9 sin 9.i) 

= R+^-riCos(9-9i), 
2R ^ 

then we can write 

gifc(K+5^-ri cos(6»-6»,)) 



Sn{R,9,t) ^ e-^^*^A 



n 

'i 



R + - Ti cos(9 - 9i] 



i=i Ji. -r 2R 

-ik(t-R) " ^2 

J2A4i)e"''-^-'^^'°<^-^^^\ (5.4) 



1=1 

which is highly accurate if R is large enough. 

For fixed t = t', a plot of the magnitude of this function will be a 
smooth closed surface in R'^ with a more or less spherical shape. This 
surface is called the wavefront at time t', and may be approximated 
locally "quite accurately" by a plane as illustrated in Figure ^ if the 
wavenumber k is not "very large". 
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Wave front 



Approximating plane 



Signal source 



Figure 2. Local approximation of wavefront by a plane. 

We now formulate our problem precisely: Is it possible, within some 
reasonable degree of accuracy, to separate Dn{a) from £)^(q;), n ^ m, 
by simply doing constant phase measurements of the signals s„, Sm in 
some small space interval far away from the signal source? 

We will use the concepts and algorithms introduced in chapter 1 to 
give an answer to this question. 



The results given in this section are very "experimental" in nature 
and were generated by implementing the Local Discriminant Basis Se- 
lection Algorithm described in Chapter 1. 

By plotting the two most discriminating coordinates in the LDB 
for a collection of signals from two different distributions, we observe 
that one class is more or less centered around the origin in a relatively 
small area whereas the other class is spread over a much larger area 
around the origin. This phenomenon is to some degree documented in 
the scatter plots shown in appendix B. This indicates that it should be 
possible to achieve a reasonable degree of discrimination between the 
two classes by using suitable "hypersphere" surfaces in constructing 
the classifier. Define the set Bn^r by 



Let 's G R"* be the point defined by the m most discriminating coor- 
dinates of some test signal s G R". Given two distributions Dx{c() ^ 
Dy{a) and m and r, we define the classifier Cm,r by 



We will concentrate on the problem of separating from Dm 
when \m — n\ — 1, since the results obtained in this case should get 
even better when \m — n\ > 1. We compute the misclassification rate 



2. Experimental results 



B, 



n,r 



{x: xl + xl-\ h < r}. 
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on both the training dataset and the test dataset for each of four dif- 
ferent real orthonormal bases: The "standard euchdean basis" (STB), 
a real Fourier-basis of type "discrete cosine basis" (DCB), the "best 
local cosine packet basis" (BLCPB) and the "best coifiet packet basis" 
(BCPB) using the coifiet of length 18. The discrete cosine basis will be 
of type IV, that is cosines evaluated at half integers in both time and 
frequency. 

We note that the time-frequency localization properties of these 
four bases are very different: The elements of the standard euclidean 
basis are perfectly localized in time but possess no localization in fre- 
quency, the elements of the discrete cosine basis are perfectly localized 
in frequency and possess no localization in time, while the elements of 
BCPB and BLCPB possess localization in both time and frequency. 

We use a training dataset of 250 signals from each class generated 
by randomly choosing points pi G Qn,i{o:) I < n < 250 in ( pTSf ), and 
a test dataset of 2500 randomly generated signals of the same type 
from each class. We set the wavenumber k to 100 and the distance of 
observation R to 10^. To not make the problem of extracting relevant 
features too easy for the LDB-algorithm, we equalize the maximum 
amplitude of the resultant signals s„ and Sn+i by setting v4„(z) = 
The real parts of the signals s„, Sn+i were sampled 2^^ = 2048 times 
with period ~ 4 ■ 10~^ units, yielding a sampling interval at the 
point of observation of length ■ 2^^ ^ 8.04 units containing about 
2^ = 128 oscillations of the signals Sn+i- The expansion depth in 
the packet bases is fixed to 8. The performance of the classifier Cm,r 
is maximized subject to the conditions: 2 < m < 10, 0.0 < r < 5.0, 
where r is an integer multiplum of 2"''. 

Given Dn, -Dn+i, we will study the performance of the LDB-algorithm 
and our particular classifier Cm,r for the cases a = 30°, a = 45°, 
a = 60°. 

2.1. The case D2(")5 -^3(')- The numbers given in boldface behind 
the error rates are the best m in the upper righthand corner and the 
best r in the lower righthand corner. The same m and r are applied to 
the test dataset, of course. The best test results for each a are written 
in bold face. 



2.2. The case D3(-), D4(-). The numbers given in boldface behind 
the error rates are the best m in the upper righthand corner and the 
best r in the lower righthand corner. The same m and r are applied to 
the test dataset, of course. The best test results for each a are written 
in bold face. 
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Method 


Error rate (%) 




Training dataset 


Test dataset 


a = 30° 


a = 45° 


a = 60° 


a = 30° 


a = 45° 


a = 60° 


STB 


35.80^27 


38.40i"39 


21.40i"43 


32.78 


33.60 


29.80 


DCB 


37.60^07 


39.60i^77 


48.60^.28 


44.52 


45.46 


49.68 


BCPB 


15.00^"55 


17.80^12 


20.20^"3i 


16.62 


22.84 


25.20 


BLCPB 


13.40^»7 


17.60^"o9 


18.80i"98 


15.94 


17.78 


19.80 



Table 1. Numerical results from testing the perfor- 



mance of Cm,r on D2{-), D^{-). 



Method 


Error rate (%) 


Cm,r on 


Training dataset 


Test dataset 


a = 30° 


a = 45° 


a = 60° 


a = 30° 


a = 45° 


a = 60° 


STB 


31.20^.8, 


31.10^1, 


lO.OOi'.oa 


11.18 


10.68 


39.91 


DCB 


38.40^.92 


45.00^.95 


40.40^^31 


38.14 


49.14 


41.68 


BCPB 


7.60^.42 


26.20^,14 


21.60^.19 


7.72 


24.94 


21.92 


BLCPB 


6.80«78 


23.80o.48 


30.80« 12 


10.56 


25.20 


31.74 



Table 2. Numerical results from testing the perfor- 



mance of Cm,r on I?3(-), D^l'). 

2.3. The case -D4(-), D^^-). The numbers given in boldface behind 
the error rates are the best m in the upper righthand corner and the 

best r in the lower righthand corner. The same m and r are applied to 
the test dataset, of course. The best test results for each a are written 
in boldface. 



Method 


Error rate (%) 




Training dataset 


Test dataset 


a = 30° 


a = 45° 


a = 60° 


q; = 30° 


q; = 45° 


a = 60° 


STB 


34.40^.63 


35.80«7e 


43.00«.63 


44.78 


44.86 


44.54 


DCB 


35.80^.61 


42.00^01 


41.20^.48 


42.32 


43.08 


45.10 


BCPB 


21.60-2 


22-80o.io 


40.80^.66 


23.62 


25.32 


43.96 


BLCPB 


15.601^42 


15.60g.32 


39.60^.73 


16.74 


18.12 


43.86 



Table 3. Numerical results from testing the perfor- 



mance of Cm,r on I?4(-), D^l'). 



2.4. The case Z}5(-), Del')- The numbers given in boldface behind 
the error rates are the best m in the upper righthand corner and the 
best r in the lower righthand corner. The same m and r are applied to 
the test dataset, of course. The best test results for each a are written 
in boldface. 
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Method 


Error rate (%) 


Cm,r on 


Training dataset 


Test dataset 


a = 30° 


a = 45° 


a = 60° 


a = 30° 


a = 45° 


a = 60° 


STB 


32.60^.65 


43.00^.58 


32.20q g2 


46.32 


46.06 


44.50 


DCB 


41.20«3o 


35.80i"74 


38.40^69 


44.74 


43.50 


44.22 


BCPB 


29.60^.45 


31.80g.6i 


41.60i% 


35.46 


33.80 


49.18 


BLCPB 


27.80»8i 


22.8O0 52 


35.20q 


32.48 


28.62 


46.72 



Table 4. Numerical results from testing the perfor- 



mance of Cm,r on D^{-), -De(-). 



2.5. The case Dio(-), -D2o(')- The numbers given in boldface be- 
hind the error rates are the best m in the upper righthand corner and 
the best r in the lower righthand corner. The same m and r are ap- 
plied to the test dataset, of course. The best test results for each a are 
written in boldface. 



Method 


Error rate (%) 




Training dataset 


Test dataset 


a = 30° 


a = 45° 


a = 60° 


a = 30° 


a = 45° 


a = 60° 


STB 


21.00U 


26.60-9 


32.00-2 


26.62 


29.22 


29.08 


DCB 


24.00-1 


24.80-4 


26.20«96 


25.90 


25.16 


26.16 


BCPB 


20.00^.67 


28.40«56 


39.60^^71 


27.56 


35.10 


34.42 


BLCPB 


22.8O1 06 


17.80-8 


20.20^.97 


24.50 


24.86 


26.72 



Table 5. Numerical results from testing the perfor- 



mance of Cm,r on -Dio(-); -D2o(-)- 



3. Conclusion 

Table |2.1| - Table show that the classifier Cm,r performs best 
on the coifiet packet and local cosine packet bases. Thus, we see that 
localization in both time/space and frequency is essential to provide the 
classifier with relevant features. We note that coifiets seem to be less 
"resistant" to overtraining than local cosines, they adapt too well to 
the training signals, that is. 



APPENDIX A 



502(1^) elements for some orthogonal FIR filters. 

Below we give the parameters 0;^, 1 < k < L, that determines 
the operators W2{ak) for 2 types of orthonormal filters of length L, 
Daubechies original shortest filters and Coifiet filters. 

1. Daubechies shortest filters. 

Daub 4 



Daub 6 



Daub 8 



ai = 0.5773502691 
a2 = -0.2679491923 



ai = 0.4122865951 
a2 = 1.831178514 
0:3 = -0.1058894200 



ai = 0.3222758836 
a2 = 1.233150027 
a3 = 3.856627874 
q;4 = -0.04600009616 



Daub 10 



ai = 0.2651451339 
a2 = 0.9398995872 
as = 2.353886784 
0:4 = 7.508378888 
as = -0.02083494630 
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Daub 12 



as 
0:4 
as 
ae 

Daub 14 

CCi 

as 
0:4 

Daub 16 

0:2 

a4 
as 
ae 

as 



0.2255061720 

0.7643296306 

1.696013010 

4.114979257 

14.28573961 

-0.009658362993 



0.1963287126 

0.6466065217 

1.333037518 

2.764759661 

7.035232916 

27.00281769 

-0.004543409641 



0.1739238836 

0.5617332940 

1.103629937 

2.074598026 

4.380557848 

12.05139151 

49.52666172 

-0.002443028170 



2. COIFLET FILTERS. 



65 



Daub 18 



Daub 20 



Coif 6 



ai = 0.1561629731 
a2 = 0.4973943657 
= 0.9452416623 
0:4 = 1.664172294 
as = 3.114016860 
ae = 6.915226655 
a-! = 20.60043019 
ag = 96.49772819 
ag = -0.001033336055 



CKi = 0.1417287200 

^2 = 0.4467987788 

as = 0.8289658876 

a4 = 1.394189716 

as = 2.402966640 

ag = 4.635603726 

a7 = 10.98508401 

ag = 35.63003753 

ag = 183.0054911 

a^o = -0.0004973444230 

2. Coiflet filters. 



ai 
as 



= -0.2152504427 
= 0.3779644639 
= -0.2152504427 



66 A. 502 (R) ELEMENTS FOR SOME ORTHOGONAL FIR FILTERS. 
Coif 12 



ai = -0.3952094767 
a2 = -0.5625481503 
eta = 0.1165449040 
ai = 1.317233974 
as = 6.198029576 
aQ = -0.04396989341 

Coif 18 

= -0.4874353702 
a2 = -1.119071133 
as = -0.2570708497 
a4 = 0.1290348165 
as = 0.4411074710 
ag = 2.215422179 
a7 = 8.338120664 
as = 15.03636438 
ag = -0.009120773147 

Coif 24 



ai = 


-0.5476023581 


0(2 — 


-1.457533881 


as = 


-0.7720754411 


a4 = 


-0.1309276144 


as = 


0.1710353887 


ae = 


0.2957793746 


ay = 


0.8070747686 


ag = 


3.126528296 


ag = 


11.27596534 


aio = 


-- 12.66598170 


an = 


-- 53.96686137 



= -0.002000409650 
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Coif 30 



r\j ^ — 
\X\ — 


U . U C/ X t: O U O t? O 


rvr, — 
'J-2 — 


— 1 7i8nnin3.'i 




— 1 iQ'inin4fiQ 


/ j — 

ct4 — 


U. UU UU UZi ±Oc/ 


rvf — 
U!5 — 


— n 1 31 fi'i39Q93 


OLq — 


n i9n'i 37301 fi 


CX-j — 


U.oD / iizDooz 


ot% = 


0.4678947012 




1.165968370 




= 4.100416655 


an = 


= 15.61099604 


Oil2 = 


= 11.59905847 


ai3 = 


= 37.56973541 


ai4 = 


= 197.1316159 




= -0.0004543371650 



A. 502 (R) ELEMENTS FOR SOME ORTHOGONAL FIR FILTERS. 



APPENDIX B 



Edge matrices for orthonormal filters 



The matrices ^Hl,- ^-re computed as follows: Each row in a ^Hl,- 
corresponds to some c or o? under the action of Eh^,. on e-?. The rows 
corresponding to c's are uniquely determined by claiming polynomials 
up to degree f — 1 mapping to polynomials of the same degree. The 
rows corresponding to (i's are not uniquely determined in this way, 
since claiming ^ — l vanishing moments on each of the rf's results in the 
matrix Ej^^ . becoming singular. This gives some "degrees of freedom" 
which were used to minimize the quantity ||E//^ .||2 + ||Ejl^^ .II2, that is 
the sum of the operator norms of the matrices, on some coarse grid. 
Remcirks: 



• The edge-matrices ^Hl,- operates on edge vectors e? resulting 
from non-normalized inner rotations. 

• Entries less than 10~^ were set to zero. 

• Superscripts, n, on the matrix entries denote multiplication by 
10"". 

• The matrix norm is the maximum column sum. 



1. Edge matrices for H^, ifg°*^ 
Applying E^d . or E^coi/ at edges leads to 

• 1 vanishing moment on do. 

• 1 vanishing moment on dn/2-i- 

• Polynomials up to and including degree 1 map to polynomials of 
the same degree. 



Grid dimension = 1, Gridsize = .1 
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B. EDGE MATRICES FOR ORTHONORMAL FILTERS 





.7719849 



9.549704 
-1.35 



.1831197 
.1047153 



1.295362 




.7865626 




-.35 
1.011213 



1.271355 




.4400402 
.9889118 





.5221843 



4.861003 
-1.95 



.7682189 
.2057189 



1.915033 




.6742919 




-.35 
1.046333 



^6 ,'■ 



1.483037 




.4960784 
.9557189 



2. Edge matrices for H[ 

Applying Yi^d ., E„coi/ leads to 



• 1 vanisliing moment on d^. 

• 2 vanisliing moments on di and dn/2-i- 

• Polynomials up to and including degree 2 map to polynomials of 
the same degree. 



Remark: The coiflet of length 8 is the one we computed numerically 
in Table |^. 

Left case: Grid dimension = 3, Gridsize = 1. 
Right case: Grid dimension = 1, Gridsize = .1 
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5.702903 



4.698626 



-7.5 



-7.876422 



-.5 



1.5 



.8137842 .7203228^ -.7748912 
.4854576 .5171221^ -.5892172 
.4590297^ 



E 



Hi,' 







4.124631 -.1195278^ 



1.186296 1.289667 -.15 
1.002116 



Hi,r 



.2635720 -.8429597 
.2424459 




-.1230332 
.2891783^ 
.9978885 



-.1888304^ -5.5 10.5 

2.666667 

-.8293807 .4881421 -1.5 



-.1067955 -.2563419 -1.203287 
-.1814515 .7167892 .4131254^ 
.3750000 



-2.519102 .7359046 
.7396251 .8283321 -.45 
1.203777 



.4445772 1.352037 .2336397 
-.3969669 .2426777 
.8307189 



3. Edge matrices for Hf^ 
Applying E^d . leads to 

• 2 vanishing moments on do and the dn/2-i- 

• 3 vanishing moments on di and the rfn/2-2- 

• Polynomials up to and including degree 3 map to polynomials of 
the same degree. 

Left case : Grid dimension = 3, Gridsize = 1. 
Right case : Grid dimension = 3, Gridsize = 1. 
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■^ff<* I 

^10'' 



( .1889731^ 
-51.61952 
.2652604^ 
\ -3.922730 

/ -.4205230^ 
-.4252506^ 
-.2725347^ 

V 



-.2458571^ 
50.72522 

-.3472944^ 
5.645330 

-.6107834^ 
-.4244118^ 







-366.9168 
.5 

.1291323^ 
-2.755913 

.2309330^ 
.1901781^ 
.1742140^ 
.2082591^ 



306.9367 ^ 

10.5 
48.01625 

1.5 / 

.5488057 \ 
.5584793 


.1007584^ / 



e: 



( 10.45706 -2.394158 

-.2181276 

20.02188 -35.08449 

\ .1417283^ 

/ .1100016 .9099426^ 

.6277529^ .9694300^ 

.1308065 

\ 



-5.5 
7.644881 
12.07371 




.5 \ 
-.2776091^ 

-.5 
1.000434 



-.7506492^ -.5847607^ \ 

-.3278640^ -.4749111^ 

.3629740^ 

.9995661 / 



4. Edge matrices for Hf^, Hlf 
Applying E^d_^ ., E^coi/ , leads to 

• 2 vanishing moments on do- 

• 3 vanishing moments on di and dn/2-\- 

• 4 vanishing moments on ^2 and dn/2-2- 

• Polynomials up to and including degree 4 map to polynomials of 
the same degree. 



Left case : Grid dimension = 6, Gridsize = 3. 
Right case : Grid dimension = 3, Gridsize = 1. 
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E 



( -9.312776 
-6.91477 
-47.71849 

-.8985454 
\ 11.10650 



16.23743 
5.863835 
24.84963 

.7639290 
-12.40799 



-10 
-3.247505 
8.319144 

-.4245146 
9.908950 



-1484.263 
-1 

.2903221 
-5.400903 



2 \ 
751.2569 

8 

103.6496 

-1 



E 



E, 



E-i 



E" 





( 


.1264939 


-.5806068^ 


-.1332401^ 


.4135265^ 


.1388293 




\ 






.2045857 


— .8716865 


.6821282 


.3779404 


.2006150 








.1144054 


— .8074226 


.2347018 


.3729639 


.1965661 








r»o/^i A -x r'S 

-.2361415 


— .0727262 


.2410276 


.4875633 


-.1600175 


4 




\ 


.5729844 





— .6966225 


.9657503 


.5300424^ 




/ 




/ 


.6614876* 


-.2352082^ 


62.56825 


-.2539657 


.2359844 


2 


\ 






— 196.6045 


298.3539 


—61.96023 


— 10.5 


.5 








—.1009016 


01 01 r\ ci a4: 

— .3181084 


.3673955 


14.35707 


-.6760348^ 






— 139.7373 


332.5374 


—210.6613 


30.02657 


-.5 






\ 


—.3170654 


.6220284 


— .3292688 


.4123728 


1.000093 


/ 




/ 


.2850594^ 


-.1403871^ 


-.3610553^ 


.1259560^ 


.1322425^ 


\ 








.2210355 


— .5899283 


— .2128227 


.8300054 


.7032465^ 










.1598260 





.2826771 





-.3750916* 
















.6965203 





.4709945* 








\ 














.9999068 


I 






( 


.3964447^ 


-2.236082 


2 


-1 










6.924950 


24.51832 


22.46071 


-97.05257 - 


-16.09946 




— 




.1434265 


— 1.288120 


—.7670091 


2 


2 








1.480799 


5.242960 


4.802957 


10.23351 


1.464548 








—.1388966 


— .4917635 


—.4504867 


—.9600 


2 ) 






( 


-.2325044 


.6560194^ 


2.279093 


.1082415^ 


-1.526436 




\ 






-.1959253 


-.3493090^ 


-.2975522 


.7647028^ 


.4346684 










.2855577 


-.5141012^ 


-.3778586 


.5405817^ 


.1133166^ 













-.6900798^ 





.2532173^ 


-.7409259^ 








\ 


-.1446773^ 





.6555962^ 


.4388504^ 


.4678566 




/ 




/ 


-6.463859 


12.73118 


6.045167 - 


-.1153122 - 


-.4242841^ 


\ 








-17.25484 


-5.83848 


24.97822 


-10.5 


.5 










1.145926 


-2.257005 


1.273759 


6.245740 - 


-.2633039^ 










-83.70978 


164.8742 


-93.04824 


9.198964 


-.5 








I 


-.4281287 


.8432407 


-.4758912 


4704796^ 


.9993761 


/ 








.8105815^ 


-.4945805^ 


-.6847096^ 


-.9052881^ 


.2054129^ 


\ 






.8381187^ 


-.2511078^ 


-.3799309^ 


-.1877050^ 


.1196993^ 








.7558504^ 





.9794265^ 


-.5689511^ 


-.2523067 


2 












.1569448 


.2142764^ 


.1113393^ 






I 











-.5104579^ 


.9980704 




/ 
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APPENDIX C 

Signal and Scatter Plots. 

1. Plots for £'2(45°), £'3(45°) 



Plot of a D2(45) signal. 




50 100 150 



Figure 1. Plot of the first 128 samples of the signals S2 
and S3. 
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Figure 2. Scatter plot of the 2 most discriminating co- 
ordinates in the standard basis for 100 signals from each 
of the classes L>2(45°) and L>3(45°). 
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Figure 3. Scatter plot of the 2 most discriminating co- 
ordinates in the discrete cosine IV basis for 100 signals 
from each of the classes -02(45°) and Ds{A5°). 
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Figure 4. Scatter plot of the 2 most discriminating co- 
ordinates in the coiflet packet basis for 100 signals from 
each of the classes .02(45°) and £)3(45°). 
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Figure 5. Scatter plot of the 2 most discriminating co- 
ordinates in the local cosine basis for 100 signals from 
each of the classes £)2(45°) and Ds{A5°). 
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2. Plots for £'3(45°), 1)4(45°) 



Plot of a D3(45) signa 
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Plot of G D4(45) signc 
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Figure 6. Plot of the first 128 samples of the signals S3 
and S4. 
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Figure 7. Scatter plot of the 2 most discriminating co- 
ordinates in the standard basis for 100 signals from each 
of the classes i:>3(45°) and ^4(45°). 
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Figure 8. Scatter plot of the 2 most discriminating co- 
ordinates in the discrete cosine IV basis for 100 signals 
from each of the classes £)3(45°) and 1)4(45°). 
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Figure 9. Scatter plot of the 2 most discriminating co- 
ordinates in the coiflet packet basis for 100 signals from 
each of the classes .03(45°) and £)4(45°). 
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Figure 10. Scatter plot of the 2 most discriminating 
coordinates in the local cosine basis for 100 signals from 
each of the classes Ds{'i5°) and £)4(45°). 
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3. Plots for L>4(45°), 1^5(45°) 
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Plot of D4(45) signal. 
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Figure 11. Plot of the first 128 samples of the signals 
S4 and S5. 
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Figure 12. Scatter plot of the 2 most discriminating 
coordinates in the standard basis for 100 signals from 
each of the classes .04(45°) and £)5(45°). 
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Figure 13. Scatter plot of the 2 most discriminating 
coordinates in the discrete cosine IV basis for 100 signals 
from each of the classes £)4(45°) and 1)5(45°). 
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Figure 14. Scatter plot of the 2 most discriminating 
coordinates in the coiflet packet basis for 100 signals from 
each of the classes .04(45°) and £)5(45°). 




-1.5 



-1.0 



-0.5 0.0 0.5 

The best LDB coordinote 



1.0 



1.5 



Figure 15. Scatter plot of the 2 most discriminating 
coordinates in the local cosine basis for 100 signals from 
each of the classes £)4(45°) and £)5(45°). 
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4. Plots for ^5(45°), L'6(45°) 




Figure 16. Plot of the first 128 samples of the signals 
S5 and sq. 
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Figure 17. Scatter plot of the 2 most discriminating 
coordinates in the standard basis for 100 signals from 
each of the classes .05(45°) and Dq{A5°). 
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Figure 18. Scatter plot of the 2 most discriminating 
coordinates in the discrete cosine IV basis for 100 signals 
from each of the classes £)5(45°) and De{^5°). 
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Figure 19. Scatter plot of the 2 most discriminating 
coordinates in the coiflet packet basis for 100 signals from 
each of the classes .05(45°) and Dq{A5°). 
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Figure 20. Scatter plot of the 2 most discriminating 
coordinates in the local cosine basis for 100 signals from 
each of the classes £)5(45°) and De{A5°). 



5. PLOTS FOR -Dio(45°), -D2o(45°) 87 

5. Plots for Dio(45°), /?2o(45°) 



Plot of D10(-15) signa 
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Figure 21. Plot of the first 128 samples of the signals 
sio and S2o- 
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Figure 22. Scatter plot of the 2 most discriminating 
coordinates in the standard basis for 100 signals from 
each of the classes £)io(45°) and £)2o(45°). 
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Figure 23. Scatter plot of the 2 most discriminating 
coordinates in the discrete cosine IV basis for 100 signals 
from each of the classes Dio{A5°) and D2o(45°). 
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Figure 24. Scatter plot of the 2 most discriminating 
coordinates in the coiflet packet basis for 100 signals from 
each of the classes £)io(45°) and £)2o(45°). 
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Figure 25. Scatter plot of the 2 most discriminating 
coordinates in the local cosine basis for 100 signals from 
each of the classes £)io(45°) and £)2o(45°). 
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APPENDIX D 



Computer Programs 

All of the numerical work discussed in this paper was carried through 
in Ansi-C or Maple. In the three sections below we give the source code 
of all important transforms and algorithms used in this paper. The code 
in the first section was, except some slight modifications, copied from 
[^. The code in the second and third section was generated by the 
author himself. 

1. Borrowed Ansi C source code 

/* Some useful functions: The discrete fourier, cosine/ sine, 
local cosine/ sine functions and their companions are copied 
(and somewhat adapted for our special use) from M.L.Wickerhauser: 
"Adapted Wavelet Analysis from Theory to Software" . */ 

10 

/* typedef complex: Define a new data structure of type "complex" 
containing two members: 

RE = Floating point number describing the 
real part of a complex number. 

IM — Floating point number describing the imaginary part 
of a complex number. */ 

typedef struct { 

double RE; 20 
double IM; 

} complex; 



/* struct interval: Define a new data structure of type "interval" , 
containing three members: 

ORIGIN = Floating point pointer to origin of an array of both 
positive and negative indexes. 

LEAST = Integer describing the least index of the array. 30 
FINAL = Integer describing the largest index of the array. */ 
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struct interval { 



double *ORIGIN; 
int LEAST; 
int FINAL; 

}; 



#define 
#define 
#define 
#define 
#define 
#define 
#define 
#define 
#define 



CCMULRE(Z1,Z2) 

CCMULIM(Z1,Z2) 

CRRMULRE(Z,YRE,YIM) 

CRRMULIM(Z,YRE,YIM) 

max(X,Y) 

min(X,Y) 

rcf(T) 

abtblock(N,L,B) 
abtlength(N,L) 



(Z1.RE*Z2.RE-Z1.IM*Z2.IM) 
(Z1.RE*Z2.IM+Z1.IM*Z2.RE) 

(Z.RE*YRE-Z.IM*YIM) 
(Z.RE*YIM+Z.IM*YRE) 



40 



((X > Y) ? X:Y) 
((X < Y) ? X:Y) 

rcfis(l,T) 
((L)*(N)+(B)*((N)»(L))) 
((N)»(L)) 



/* makeinterval: Allocate an interval data structure 
and assign its data array */ 



50 



struct interval makeinterval(double *DATA, int LEAST, int FINAL) 
{ 

int LENGTH, K; 

struct interval SEG; 

LENGTH = 1+FINAL-LEAST; 
if (LENGTH > 0) { 
SEG.ORIGIN = (double *)calloc(LENGTH, sizeof (double)); 60 
SEG.ORIGIN -= LEAST; 
if (DATA != NULL) { 
for (K = LEAST; K <= FINAL; K++) { 
SEG.ORIGIN[K] = DATA[K-LEAST]; 

} 

} 

} 

SEG.LEAST = LEAST; 
SEG.FINAL = FINAL; 

return SEG; 70 



} 



/* freeinterval: Deallocate an interval data structure 
and its data array */ 



double *£reeinterval(struct interval *SEG) 
{ 

if (SEG != NULL) { 



if ((SEG->ORIGIN) != NULL) { 



(SEG->ORIGIN) += SEG->LEAST; 
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free(SEG->ORIGIN); 

} 

free(SEG); 

} 

return NULL; 



/* br: Return the input integer bit-reversed */ 

int br(int N, int L0G2LEN) 
{ 

int U, J; 
U = N.S^1; 

for (J = 1; J < L0G2LEN; J++) { 
N »= 1; 

U <<= 1; 
U += N&il; 

} 

return U; 

} 

/* bitrevd: Permute to a disjoint array by bit-reversing 
the indices */ 

void bitrevd(coniplex *OUT, complex *IN, int Q) 
{ 

int M, N, U; 
M = 1«Q; 

for (N 0; N < M; N++) { 
U = br(N, Q); 
OUT[U] = IN[N]; 

} 

} 

/* bitrevi: Permute an array in place via index bit-reversal */ 

void bitrevi (complex *X, int Q) 
{ 

int M, N, U; 

complex TEMP; 

M = 1«Q; 

for (N = 1; N < (M-1); N++) { 
U = br(N, Q); 
if (U > N) { 
TEMP = X[N]; 
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X[N] = X[U]; 130 
X[U] = TEMP; 

} 

} 

} 

/ * fflomega: Compute table of sines and cosines for DFT * / 

void fFtomega(complex *W, int M) 
{ 

int K; 140 

double FACTOR; 

FACTOR = (double) (-PI/M); 
if (M < 0) { 

M = -M; 

} 

for (K = 0; K < M; K++) { 
W[K].RE = cos(K*FACTOR); 

W[K].IM = sin(K*FACTOR); 150 

} 

} 

/ * fftproduct: Product of sparse matrices for DFT * / 

void fftproduct(complex *F, int Q, complex *W) 

{ 

int K, N, J, M, Nl, B; 

complex TMP; 160 

N = 1«Q; 
K = Q; 

while (K > 0) { 

K -= 1; 

Nl = N>>K; /* Block size */ 
M = Nl/2; /* Butterfly size */ 
B = 0; 

while (B < N) { 

TMP.RE = F[B+M].RE; 170 

TMP.IM = F[B+M].IM; 

F[B+M].RE = F[B].RE - TMP.RE; 

F[B+M].IM = F[B].IM - TMP.IM; 

F[B].RE += TMP.RE; 

F[B].IM += TMP.IM; 



for (J = 1; J < M; J++) { 
TMP.RE = CCMULRE(F[B+M+J], W[J*(N/N1)]); 



1. BORROWED ANSI C SOURCE CODE 

TMP.IM = CCMULIM(F[B+M+J], W[J*(N/N1)]); 
F[B+M+J].RE = F[B+J].RE - TMP.RE; 
F[B+M+J].IM = F[B+J].IM - TMP.IM; 
F[B+J].RE += TMP.RE; 
F[B+J].IM += TMP.IM; 

} 

B += Nl; 

} 

} 

} 

/* dct4omega: Compute table of sines and cosines for 
DCT-IV, DST-IV*/ 

void dct4oinega(double *COS, double *SIN, int M) 
{ 

int K; 

double FACTOR; 

FACTOR = PI/(2.0*M); 

for (K = 0; K < M; K++) { 

COS[K] = cos(K*FACTOR); 

SIN[K] = sin(K*FACTOR); 

} 

} 

/* dctnorm,al: Normalization for unitary L-point DCT 
and DST (T is log L) */ 

void dctnormal(double *Z, int L, int T) 

{ 

double NORM; 
int K; 

if ((T%2) != 0) 

NORM = 0.5/(l«((T+l)/2)) ; 
else 

NORM = 0.5/sqrt(l«(T+l)); 
for (K = 0; K < L; K++) { 
Z[K] *= NORM; 

} 

} 

/* dct4-: (in place) Unitary DCT-IV. 
DBS! Q is log of vectorlength * / 



void dct4(double *X, int Q) 
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{ 

int N, K; 

complex *F, *W, U, TMP; 
double *C, *S, *Y; 
N = 1«Q; 

F = (complex *)calloc(2*N, sizeof (complex)); 
W = (complex *)calloc(N, sizeof (complex)); 
C = (double *)calloc(N, sizeof (double)); 
S = (double *)calloc(N, sizeof (double)); 

F = (double *)calloc(N, sizeof (double)); */ 

dct4omega(C, S, N); 
F[0].RE = X[0]; 
F[N].IM = X[N-1]; 

for (K = 1; K < N; K++) { 
F[K].RE = X[K] * C[K]; 
F[K].IM = -X[K] * S[K]; 
F[2*N-K].RE = X[K-1]*C[K]; 
F[2*N-K].IM = X[K-1]*S[K]; 

} 

bitrevi(F, Q+1); 
fftomega(W, N); 

fftproduct(F, Q+1, W); 
U.RE = cos(-PI/(4.0*N)); 
U.IM = sin(-PI/(4.0*N)); 

TMP.RE = F[0].RE + CRRMULRE(F[2*N-1], C[l], S[l]); 
TMP.IM = F[0].IM + CRRMULIM(F[2*N-1], C[l], S[l]); 
X[0] = CCMULRE(TMP, U); 

TMP.RE = CRRMULRE(F[N-1], C[N-1], -S[N-1]) - F[N].IM; 
TMP.IM = CRRMULIM(F[N-1], C[N-1], -S[N-1]) + F[N].RE; 
X[N-1] = CCMULRE(TMP, U); 

for (K = 1; K < (N-1); K++) { 

TMP.RE = CRRMULRE(F[K], C[K], -S[K]) + 
CRRMULRE(F[2*N-K-1], C[K+1], S[K+1]); 

TMP.IM = CRRMULIM(F[K], C[K], -S[K]) + 
CRRMULIM(F[2*N-K-1], C[K+1], S[K+1]); 

X[K] = CCMULRE(TMP, U); 

} 

free(W); 
free(C); 
free(S); 
free(F); 

dctnormal(X, N, Q); 
/* return(Y); */ 
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/* rcfis: Iterated sine rising cutoff function */ 280 

double rcfis(int N, double T) 
{ 

int I; 

if (T > -1.0) { 
if (T < 1.0) { 
for (I = 0; I < N; I++) { 
T = sin(0.5*PI*T); 

} 290 
T = sin(0.25*PI*(1.0 + T)); 

} 

else 
T = 1.0; 

} 

else 

T = 0.0; 
return T; 

} 

300 

/* rcfmidp: Rising cutoff function sampled between gridpoints */ 



void rcfniidp(struct interval R) 
{ 

int J; 

double X, DX; 

X = 0.5/(R.FINAL + 1.0); 

DX = 1.0/(R.FINAL + 1.0); 

for (J = 0; J <= R.FINAL; J++) { 

R.ORIGIN[J] = rcf(X); 

R.0RIGIN[-J-1] = rcf(-X); 

X += DX; 

} 



310 



} 

/* fdcn: (midpoint) Fold disjoint cosine negative */ 

void fdcn(double *ONEG, int STEP, double *INEG, 320 
double *IPOS, int N, struct interval RISE) 

{ 

int K; 



for (K = -N; K <= (RISE.LEAST - 1); K++) { 
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*(ONEG+K*STEP) = *(INEG+K); 

} 

for (K = RISE. LEAST; K -1; K++) { 
*(ONEG+K*STEP) = RISE.0RIGIN[-1-K]*(*(INEG+K)) - 
RISE.0RIGIN[K]*(*(IP0S-1-K)); 330 

} 

} 



/* fdcp: (midpoint) Fold disjoint cosine positive */ 

void fdcp(double *OPOS, int STEP, double *INEG, 
double *IPOS, int N, struct interval RISE) 

{ 

int K; 340 

for (K = 0; K RISE.FINAL; K++) { 
*(OPOS+K*STEP) = RISE.ORIGIN[K]*(*(IPOS+K)) + 
RISE.0RIGIN[-1-K]*(*(INEG-1-K)); 

} 

for (K = (RISE.FINAL + 1); K (N-1); K++) { 
*(OPOS+K*STEP) = *(IPOS+K); 

} 

} 

350 

/* fdsn: (midpoint) Fold disjoint sine negative */ 

void fdsn(double *ONEG, int STEP, double *INEG, 
double *IPOS, int N, struct interval RISE) 

{ 

int K; 

for (K = -N; K <= (RISE.LEAST - 1); K++) { 

*(ONEG+K*STEP) = *(INEG+K); 
} 360 
for (K = RISE.LEAST; K <= -1; K++) { 

*(ONEG+K*STEP) = RISE.0RIGIN[-1-K]*(*(INEG+K)) + 
RISE.0RIGIN[K]*(*(IP0S-1-K)); 

} 

} 

/* fdsp: (midpoint) Fold disjoint sine positive */ 

void fdsp(double *OPOS, int STEP, double *INEG, 

double *IPOS, int N, struct interval RISE) 370 

{ 

int K; 



for (K = 0; K <= RISE.FINAL; K++) { 
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*(OPOS+K*STEP) = RISE.ORIGIN[K]*(*(IPOS+K)) - 
RISE.0RIGIN[-1-K]*(*(INEG-1-K)); 

} 

for (K = (RISE.FINAL + 1); K <= (N-1); K++) { 
*(OPOS+K*STEP) = *(IPOS+K); 

} 

} 

/* localcosine: (in place) (adapted dyadic) local cosine analysis 
(with) fixed folding applying a fixed cutoff function "RISE" 
sampled between gridpoints resulting in L+1 different 
representations of the signal. 

We have: 

N = The length of the signal to be analysed. 

L = The depth of expansion in the dictionary of bases. 

PARENT = Pointer to array of length N*(L+1) containing 
the signal to be analysed in its first N locations. */ 

void localcosine(double *PARENT, int N, int L) 

{ 

int NP, NC, LEVEL, PBLOCK, RADIUS; 

double *MIDP, *CHILD; 

struct interval RISE; 

RADIUS = (N»(L+1))-1; 
if (RADIUS <= 0) 
printf(" Negative folding radius in lcaff\n"); 

RISE = makeinterval(NULL, -RADIUS, RADIUS-1); 

rcfniidp(RISE); 

NP = N; 

for (LEVEL = 0; LEVEL < L; LEVEL++) { 
NC = NP/2; 

for (PBLOCK = 1; PBLOCK <= (1«LEVEL); PBLOCK++) { 
MIDP = PARENT + NC; 
CHILD = MIDP + N; 
fdcii(CHILD, 1, MIDP, MIDP, NC, RISE); 
fdcp(CHILD, 1, MIDP, MIDP, NC, RISE); 
dct4(PARENT, intlog2(NP)); 
PARENT += NP; 

} 
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NP = NC; 

} 

for (PBLOCK = 1; PBLOCK <= (1«L); PBLOCK++) { 
dct4(PARENT, intlog2(NP)); 
PARENT += NP; 

} 

freeinterval(,SiRISE); 430 



} 
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101 



#include <stdio.h> 
#include <math.h> 
#include <stdlib.h> 
#include <time.h> 
#include <assert.h> 



/* basisnode: Define a new data structure of type basisnode 
containing five members: 

LEVEL = The level coordinate of a subspace in the dictionary tree. 10 
BLOCK = The position of a subspace counted from, the left. 
POS = The coordinate of a basisfunction in the subspace. 
COST = The cost/ discrimination number. 
TAG = YES or NO. */ 



struct basisnode { 
int LEVEL; 
int BLOCK; 
int POS; 

double COST; 20 
int TAG; 

}; 



#define 


YES 


1 


#define 


NO 





#define 


PI 


3.1415926535 


#define 


SQ2 


1.4142135623 


#define 


SQH 


0.7071067811 


#define 


MAXINT 


2147483647 


#define 


DIVISOR 


127773 


#define 


HIGH 


2836 


#define 


LOW 


16807 


#define 


SEEDFACTOR 


76928675 


#define 


CLASS_1 


10 


#define 


CLASS_2 


20 


#define 


K_l 


pow(10. 0,2.0) 


#define 


K_2 


pow(10.0,2.0) 


#define 


A_l 


0.100 


#define 


A_2 


0.050 


#define 


RDIST 


10000 


#define 


expand 


localcosine 


#define 


s_measure(IN_l,IN_2,N) 


lp_norm_p(IN_l,IN_2,N,2.0) 


#define 


makcsignal (IN,N,K,A,M) 


make_samples_of_cosine_signal( 


IN,N,K,A,RDIST,16,M,1.0/8.0,PI/4) 





/* Allocate a basisnode data structure an assign its content */ 
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struct basisnode *makebasisnode(int LEVEL, int BLOCK, 

int POS, double COST, int TAG) 

{ 

struct basisnode *BASISNODE; 

BASISNODE = (struct basisnode *)malloc(sizeof(struct basisnode)); 

BASISNODE->LEVEL = LEVEL; 
BASISNODE->BLOCK = BLOCK; 
BASISNODE->POS = POS; 
BASISNODE->COST = COST; 
BASISNODE->TAG = TAG; 

return(BASISNODE); 

} 

/* Deallocate an array of length N of pointers 
to basisnode data structures */ 

void freebasisnodes(struct basisnode **A, int N) 
{ 

int K; 

for (K = 0; K < N; K++) { 
free(A[K]); 

} 

free(A); 

} 

/* ladder: (disjoint) filter a real array of length N into 

lowpass and highpass- coefficients by the "rotation" algorithm 
applying an orthonormal filter of length 2*L, periodizing 
at endpoints. 

We have: 

IN = Pointer to array of numbers to be filtered. 

OUT = Pointer to array that results from filtering IN. 

ALPHA = Pointer to array of "rotation" coefficients including 
a normalization factor. 

N = The length of IN = The length of OUT. 
L = Half the length of the filter. */ 
void ladder(double *OUT, double *IN, 
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double *ALPHA, int N, int L) 

{ 

int LEVEL, TRANSLATION, K; 
double *WORK; 

WORK = (double *)calloc(N, sizeof (double)); 
for (K = 0; K < N; K++) 
WORK[K] = IN[K]; 

for (LEVEL = 0; LEVEL < L; LEVEL++) { 
TRANSLATION = (LEVEL%2); 

for (K = 0; K < N/2; K++) { 

*(OUT+TRANSLATION + 2*K) = 
WORK[TRANSLATION + 2*K] - ALPHA[LEVEL]* 
WORK[(TRANSLATION + 1 + 2*K)%N]; 

*(OUT+(TRANSLATION + 1 + 2*K)%N) = 
ALPHA[LEVEL]*WORK[TRANSLATION + 2*K] + 
WORK[(TRANSLATION + 1 + 2*K)%N]; 

WORK[TRANSLATION + 2*K] = 

*(OUT+TRANSLATION + 2*K); 
WORK [(TRANSLATION + 1 + 2*K)%N] = 

*(OUT+(TRANSLATION + 1 + 2*K)%N); 

} 

} 

for (K = 0; K < N/2; K++) { 

/* Store the lowpass and highpass coefficients in 
separate blocks. */ 

if (L%2 == 0) { 
*(OUT+K) = W0RK[2*K +1]*ALPHA[L]; /* c's */ 

*(0UT+K+N/2) = W0RK[2*K]*ALPHA[L]; /* d's */ 

} 

else { 

*(OUT+K) = W0RK[2*K]*ALPHA[L]; /* c's */ 

*(0UT+K+N/2) = W0RK[2*K+1]*ALPHA[L]; /* d's */ 

} 

} 

free(WORK); 



/* waveletpacket: (in place) waveletpacket analysis resulting in L+1 
different representations of the signal. 
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We have: 

N = The length of the signal to be analysed. 

L = The depth of expansion in the dictionary of bases. 

PARENT = Pointer to array of length N*(L+1) containing 
the sampled lowpass coefficients of the signal to be analysed 
in its first N locations. 

ALPHA = Pointer to array containing in the following order: 

1. Half the length of the orthonormal filter. 

2. The rotation coefficients corresponding to the filter. 

3. A scaling factor. */ 

void waveletpacket (double *PARENT, int N, int L) 
{ 

FILE *FILTER; 

int LEVEL, BLOCK, FILEN, NP, M, K; 
double * CHILD, * ALPHA; 

if(( FILTER = fopen("/users/stud/eirikf /programs/hovedf ag/ 
f ilters/laddercoif 18","r")) 
== NULL) { 

printf("ERROR in opening FILTERf ile\n"); 
exit(l); 

} 

fscanf (FILTER, "7.d", .SdFILEN); 

ALPHA = (double *)calloc(FILEN+l, sizeof(double)); 
for (K = 0; K <= FILEN; K++) { 
fscanf (FILTER, "%lf", &ALPHA[K]); 

} 

fclose(FILTER); 

for (LEVEL = 0; LEVEL < L; LEVEL++) { 
NP = N»LEVEL; 

for (BLOCK = 0; BLOCK < (1«LEVEL); BLOCK++) { 
CHILD = PARENT+N; 

laddcr(CHILD, PARENT, ALPHA, NP, FILEN); 
PARENT += NP; 

} 

} 
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£ree( ALPHA); 

} 

/* my_qsort: sort v[left]. . .v[right] into decreasing order. 

(This function is copied from Kernighan/ Ritchie: 200 
"The Ansi C programming language" second edition p. 120.) */ 

void my_qsort(void *v[], int left, int right, 
int (*comp)(void *, void *)) 

{ 

int i, last; 

void my_swap(void *v[], int, int); 

if (left >= right) 210 
return; 

my_swap(v, left, (left + right)/2); 

last = left; 

for (i = lcft+1; i <= right; i++) { 
if ((*comp)(v[i], v[left]) < 0) 
my_swap(v, ++last, i); 

} 

my_swap(v, left, last); 

my_qsort(v, left, last— 1, comp); 220 
my_qsort(v, last+1, right, comp); 

} 

/* myswap: Exchange pointers within an array */ 

void my_swap(void *v[], int i, int j) 
{ 

void *temp; 

temp = v[i]; 230 
v[i] = vp]; 
vp] = temp; 

} 

/* comp_basisnode-by-Cost: Compare the discrimination measures 
of two basisnode data structures */ 

int comp_basisnode_by_cost(struct basisnode *X, 

struct basisnode *Y) 
{ 240 
if ((X->COST) > (Y->COST)) 

return(— 1); 
else if ((X->COST) < (Y->COST)) 
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return(l); 
else 
return(O); 

} 

/* my_rand: Generate a pseudo random floating point 

number between and 1. */ 250 

double my_rand(void) 

int high, low, temp; 
static int SEED; 

if (SEED < 1 I I SEED > MAXINT) 
SEED = ( (int) time(NULL)*SEEDFACTOR ) % MAXINT; 

260 

high = SEED/DIVISOR; 

low = SEED % DIVISOR; 

temp = LOW*low - HIGH*high; 

SEED = (temp > ? temp : temp + MAXINT); 

return( (double)SEED /MAXINT) ; 

} 

/* lp_norm_p: Compute the (l~P difference) ~P of two arrays 
of length N */ 

270 

double lp_norm_p(double *IN_1, double *IN_2, int N, double P) 

{ 

double T, X; 
int K; 

T = 0; 

if (IN_2 == NULL) { 
for (K = 0; K < N; K++) { 
X = fabs(*(IN_l+K)); 280 
if (X > 0) 
T +=exp(P*log(X)); 

} 

} 



else { 

for (K = 0; K < N; K++) { 
X = fabs(*(IN_l+K)-(*(IN_2+K))); 
if ( X > 0) 
T += exp(P*log(X)); 

} 

} 



290 
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return T; 

} 

/* modf-2PI: Compute remainder of input modulo 2*PI. 
(This is because of the bad accuracy of the 
standard math trigonometric functions for large aruments) */ 

double modf_2PI(double X) 300 

{ 

int K; 

double T, S; 

T = X; 

S = 2.0*PI; 
if (T > 0) { 

while (T > 0) 310 
T -= S; 

T += S; 

} 

else { 
while (T < 0) 

T += S; 

} 

return(T); 

320 

} 

/* set_tag_to_zero: Set the tags of all children nodes of a 
basisnode data type, to NO. */ 

void set_tag_to_zero(struct basisnode **BASIS, int L, 
int LEVEL, int BLOCK) 

{ 

int PARENT, LEFTCHILD, RIGHTCHILD; 

330 

PARENT = (1«LEVEL) - 1 + BLOCK; 
LEFTCHILD = 2*PARENT + 1; 
RIGHTCHILD = LEFTCHILD + 1; 
BASIS [LEFTCHILD] -> TAG = NO; 
BASIS [RIGHTCHILD] -> TAG = NO; 
if (LEVEL < (L - 1)) { 

set_tag_to_zero(BASIS, L, LEVEL+1, 2*BL0CK); 

set_tag_to_zero(BASIS, L, LEVEL+1, 2*BL0CK + 1); 

} 

340 

} 
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/* make_samples_of_cosines: Returns an array of samples of 
a specific (cosine)function. 



We have: 



SIGNAL = Array of length N*L 



N 



The length of the sampling. 



350 



K 



The wavenumber of the transmitted wave(s). 



A 



Amplitude (factor) of signal. 



R 



The distance from origo to the interval arc of sampling. 



P 



The number of samples pr. wavelength. 



M 



The number of signal transmitters 



360 



T = Starting point of sampling. 

W = Width of sector containing the signal emitter 
as a fraction of 2.0PI. */ 

void make_samples_of_cosine_signal(double *SIGNAL, int N, double K, 

double A, double R, int P, int M, 
double W, double T) 
{ 370 
int I, J, L; 

double RADIUS, ANGLE, SAMPL, AMP, ARG, X, Y, *WORK; 
SAMPL = 2.0*PI/(K*P); 

WORK = (double *)calloc(N, sizeof (double)); 
RADIUS = 0; 

L = 1; 380 

for (J = 1; J <= M; J++) { 
while (RADIUS < 1.0) 
RADIUS = 10.0*my_rand(); 

ANGLE = 2.0*PI*(L*1.0/M + my_rand()*W); 

AMP = A; 

X = K*RADIUS*RADIUS/(2.0*R); 
Y = K*RADIUS; 



390 
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for (I = 0; I < N; I++) { 
ARG = T + I*SAMPL; 

/* WORKfl] += AMP*(1.0I R)*cos(modf-2PI(X - 

Y*cos(ARG - ANGLE))); */ 
WORKp] += AMP*cos(modf_2PI(X - Y*cos(ARG - ANGLE))); 
SIGNALp] = WORKp]; 

} 

L += 1; 
RADIUS = 0; 

} 400 
free(WORK); 

} 

/* time-frequency-energy-map: Constructs the training 

signals and their expansions into some dictionary of orthonormal 

wavelet bases or local trigonometric bases. The resulting 
"time-frequency energy map" for each class, and the basis 
coordinates of each training signal are stored in files. 

410 

We have: 

N = The length of the training signals. (Power of 2) 

L = The depth of expansion in the dictionary tree of bases. 

NS = The number of training signals in each class. 

makesignal = Function that constructs a signal 

by some sampling procedure. 420 

expand = Function that expands the signal into 

a dictionary of the desired orthonormal basis functions. */ 

void tinie_frequency_energy_inap(int N, int L, int NS, FILE *TF1, 

FILE *TF2, FILE *SIG1, FILE *SIG2) 

{ 

double L0CENERGY_1, L0CENERGY_2, CLASSN0RMP_1, 
CLASSN0RMP_2,*W0RK_1, *W0RK_2, 

*GAMMA_1, *GAMMA_2; 430 
int LEVEL, BLOCK, POS, K, J; 

GAMMA_1 = (double *)calloc(N*(L+l), sizeof(double)); 
GAMMA_2 = (double *)calloc(N*(L+l), sizeof(double)); 
W0RK_1 = (double *)calloc(N*(L+l), sizeof (double)); 
W0RK_2 = (double *)calloc(N*(L+l), sizeof (double)); 
L0CENERGY_1 = 0.0; 
L0CENERGY_2 = 0.0; 
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CLASSN0RMP_1 = 0.0; 
CLASSN0RMP_2 = 0.0; 
fprintf(SIGl, "7.d\n", N); 
fprmtf(SIGl, "y.d\n", L); 
fprintf(SIGl, "%d\n", NS); 
fprintf(SIG2, "7.d\n", N); 
fprintf(SIG2, "7.d\n", L); 
fprintf(SIG2, "y.d\n", NS); 

for (J = 1; J <= NS; J++) { 
makesignal(WORK_l, N, K_l, A_l, CLASS_1); 
makesignal(W0RK_2, N, K_2, A_2, CLASS_2); 

/* make-samples-of-l(WORK_l, N*(L+1)); 

make-samples-of-Oscillating-fcn(WORK_2, N*(L+1), 32); */ 

CLASSN0RMP_1 += lp_norin_p(WORK_l, NULL, N, 2.0); 
CLASSN0RMP_2 += lp_norni_p(WORK_2, NULL, N, 2.0); 
expand(WORK_l, N, L); 
expand(W0RK_2, N, L); 

for (K = 0; K < N*(L+1); K++) { 
fprintf(SIGl, "y.leXn", W0RK_1[K]); 
fprintf(SIG2, "/.leXn", W0RK_2[K]); 

} 

for (LEVEL = 0; LEVEL L; LEVEL++) { 
for (BLOCK = 0; BLOCK < (1«LEVEL); BLOCK++) { 
for (POS = 0; POS < abtlength(N, LEVEL); POS++) { 

L0CENERGY_1 = 
(*(W0RK_1 + abtblock(N,LEVEL,BLOCK) + POS))* 
(*(W0RK_1 + abtblock(N,LEVEL,BLOCK) + POS)); 

L0CENERGY_2 = 
(*(W0RK_2 + abtblock(N,LEVEL,BLOCK) + POS))* 
(*(W0RK_2 + abtblock(N,LEVEL,BLOCK) + POS)); 

*(GAMMA_1 + abtblock(N,LEVEL,BLOCK) + POS) 

+= L0CENERGY_1; 
*(GAMMA_2 + abtblock(N,LEVEL,BLOCK) + POS) 

+= L0CENERGY_2; 

if (J NS) { 
*(GAMMA_1 + abtblock(N,LEVEL,BLOCK) + POS) 

/= CLASSN0RMP_1; 
*(GAMMA_2 + abtblock(N,LEVEL,BLOCK) + POS) 

/= CLASSN0RMP_2; 
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} 

} 

} 



490 



} 

} 

fprintf(TFl, 
fprintf(TFl, 
fprintf(TF2, 
fprintf(TF2, 




N); 

L); 

N); 

L); 



for (K = 0; K < N*(L+1); K++) { 
fprintf(TFl, "Ue\n" , GAMMA_1[K]); 
fprintf(TF2, "y.le\n", GAMMA_2[K]); 



500 



} 



frcc(WORK_l); 
free(W0RK_2); 
£ree(GAMMA_l); 

frcc(GAMMA_2); 
fclosc(TFl); 

fclose(TF2); 510 

fclose(SIGl); 

fclose(SIG2); 

} 

/* Idb-search-and-sori: Do a LDB search and sorting procedure 
on two given binary tree arrays X and Y of time-frequency 
energy maps resulting from a collection of wavelet packet 
dictionarys or local cosine/ sine dictionarys of training 
signals belonging to two distinct classes, maximizing some 

given discriminant measure on the time-frequency energy maps. 520 
We have: 

TF_1 = Pointer to the file of time-frequency energy 
maps resulting from the set of expanded training signals 
belonging to a certain class, arranged level by level, 
block by block. 

TF_2 = Pointer to the file of time-frequency energy 

maps resulting from the set of expanded training signals 530 
belonging to some other class, arranged level by level, 
block by block. 

OUT = Pointer to file of output data. 



N = The length of each training signal. 
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L = The depth of expansion of each signal in 
the dictionarys. */ 

540 

void ldb_search_and_sort(FILE *TF_1, FILE *TF_2, FILE *OUT) 
{ 

int LEVEL, BLOCK, BLOCKLEN, POS, PARENT, 
NODE, CHILDNODE, N, L, K; 

double DELTA, DELTACHI, *W0RK_1, *W0RK_2; 

struct basisnode **BASIS, **TEMP; 

550 

fscanf(TF_l, "7.d", i^N); 
fscanf(TF_l, "7.d", &L); 

fscanf(TF_2, "%d", &N); /* Advancement of FILE pointer */ 
fscanf(TF_2, "%d", &L); /* Advancement of FILE pointer */ 

W0RK_1 = (double *)calloc(N*(L+l), sizeof(double)); 
W0RK_2 = (double *)calloc(N*(L+l), sizeof (double)); 
TEMP = (struct basisnode **)calloc(l<<(L+l), 

sizeof (struct basisnode *)); 
BASIS = (struct basisnode **)calloc(N, 560 

sizeof(struct basisnode *)); 

for (K = 0; K < N*(L+1); K++) { 
fscanf(TF_l, "7.1f ", &W0RK_1[K]); 
fscanf(TF_2, ""/,lf", &W0RK_2[K]); 

} 

BLOCKLEN = abtlength(N,L); 

for (BLOCK = 0; BLOCK < (1«L); BLOCK++) { 
NODE = (1«L) - 1 + BLOCK; 570 
DELTA = s_nieasure(WORK_l + L*N + BLOCK*BLOCKLEN, 

W0RK_2 + L*N + BLOCK*BLOCKLEN, 

BLOCKLEN); 

TEMP[NODE] = makebasisnode(L, BLOCK, 0, DELTA, YES); 

} 

for (LEVEL = L-1; LEVEL >= 0; LEVEL ) { 

BLOCKLEN = abtlength(N,LEVEL); 

for (BLOCK = 0; BLOCK < (1«LEVEL); BLOCK++) { 580 
PARENT = abtblock(N, LEVEL, BLOCK); 
NODE = (1<<LEVEL) - 1 + BLOCK; 
CHILDNODE = 2*N0DE + 1; 
DELTA = s_measure(WORK_l+PARENT, 

W0RK_2+PARENT, BLOCKLEN); 
DELTACHI = TEMP [CHILDNODE] -> COST + 
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TEMP[CHILDNODE + 1]->C0ST; 
if (DELTA >= DELTACHI) { 
set_tag_to_zero(TEMP, L, LEVEL, BLOCK); 
TEMP[NODE] = makebasisnode(LEVEL, BLOCK, 0, 

DELTA, YES); 

} 

else 

TEMP[NODE] = makebasisnode(LEVEL, BLOCK, 0, 

DELTACHI, NO); 

} 

} 

K = 0; 

for (NODE = 0; NODE < ((1«(L+1))-1); NODE++) { 
if (K < N) { 
if (TEMP[NODE]->TAG == YES) { 
LEVEL = TEMP [NODE] ->LEVEL; 
BLOCK = TEMP[NODE]->BLOCK; 
for (POS = 0; POS < (N»LEVEL); POS++) { 
DELTA = s_nieasurc(WORK_l + 

abtblock(N, LEVEL, BLOCK) + POS, 
W0RK_2 + 

abtblock(N, LEVEL, BLOCK) + POS, 1); 
BASIS[K] = niakebasisnode(LEVEL, BLOCK, POS, 
DELTA, YES); 

K += 1; 

} 

} 

} 

} 

my_qsort((void **)BASIS, 0, N-1, 

(int (*)(void *, void *))comp_basisnode_by_cost); 
fprintf(OUT, "7.d\ny.d\n" , N, L); 

for (K = 0; K < N; K++) 
fprintf(OUT, "%d %d %d y.le\n\ii", 

BASIS[K]->LEVEL, BASIS[K]->BLOCK, 
BASIS[K]->POS, BASIS[K]->COST); 

free(WORK_l); 
£ree(W0RK_2); 

frccbasisnodcs(TEMP, 1<<(L+1)); 

freebasisnodes(BASIS, N); 

fclose(TF_l); 

fclosc(TF_2); 

fclose(OUT); 
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/* make2bestcoordinates: Make a set of signals from each class 
and write out the coordinates of the two most discriminating 
LDB basis elements. 

We have: 640 

NS = The number of signals to be generated. 

LDB = Pointer to file containing the LDB data. 

CLICI = Pointer to file containing the coordinates of the 
best LDB basis element of CLASS 1. 

CL1C2 = Pointer to file containing the coordinates of the 

second best LDB basis element of CLASS 1. 650 

CL2C1 = Pointer to file containing the coordinates of the 
best LDB basis element of CLASS 2. 

CL2C2 = Pointer to file containing the coordinates of the 
second best LDB basis element of CLASS 2. */ 

void make2bestcoordinates(int NS, FILE *LDB, FILE *CL1C1, 

FILE *CL1C2, FILE *CL2C1, FILE *CL2C2) 
{ 660 
double *SIGNAL_1, *SIGNAL_2, XCOORDCLl, YCOORDCLl, 
XC00RDCL2, YC00RDCL2, NORMPl, N0RMP2, DISC; 

struct basisnode **BASIS; 

int N, L, LEVEL, BLOCK, POS, K, J; 

fscanf(LDB, "%d", &N); 
fscanf(LDB, "%d", &L); 

BASIS = (struct basisnode **)calloc(2, sizeof (struct basisnode *)); 670 
SIGNAL_1 = (double *)calloc(N*(L+l), sizeof(double)); 
SIGNAL_2 = (double *)calloc(N*(L+l), sizeof (double)); 

for (K = 0; K < 2; K++) { 
if (K < N) { 
fscanf(LDB, "7.d", &LEVEL); 
fscanf(LDB, "7.d", &BLOCK); 
fscanf(LDB, "%d", .S^POS); 
fscanf(LDB, "7.1f ", &DISC); 

BASIS[K] = niakebasisnode(LEVEL, BLOCK, POS, DISC, 1); 680 

} 

} 

fprintf(CLlCl, "7.d\n", NS); 
fprintf(CLlC2, "7.d\n", NS); 
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fprintf(CL2Cl, "'/.d\n", NS); 
fprintf(CL2C2, "%d\ii", NS); 

for (J = 0; J < NS; J++) { 
niakesignal(SIGNAL_l, N, K_l, A_l, CLASS_1); 
makesignal(SIGNAL_2, N, K_2, A_2, CLASS_2); 
NORMPl = lp_norni_p(SIGNAL_l, NULL, N, 2.0); 
N0RMP2 = lp_norni_p(SIGNAL_2, NULL, N, 2.0); 
expand(SIGNAL_l, N, L); 
expand(SIGNAL_2, N, L); 

XCOORDCLl = *(SIGNAL_1 + 

abtblock(N, BASIS[0]->LEVEL, BASIS [0]->BLOCK) + 

BASIS[0]->POS); 
XC00RDCL2 = *(SIGNAL_2 + 

abtblock(N, BASIS[0]->LEVEL, BASIS [0]->BLOCK) + 

BASIS[0]->POS); 
YCOORDCLl = *(SIGNAL_1 + 

abtblock(N, BASIS[1]->LEVEL, BASIS[1]->BL0CK) + 

BASIS[1]->P0S); 
YC00RDCL2 = *(SIGNAL_2 + 

abtblock(N, BASIS[1]->LEVEL, BASIS[1]->BL0CK) + 

BASIS[1]->P0S); 
fprintf(CLlCl, "7.1f\n", XCOORDCLl); 
fprintf(CLlC2, "7.1f\n", YCOORDCLl); 
fprintf(CL2Cl, "7,lf\n", XC00RDGL2); 
fpriiitf(CL2C2, "•/,lf\n", YG00RDCL2); 

} 

freebasisnodes(BASIS, 2); 

£ree(SIGNAL_l); 

£ree(SIGNAL_2); 

fclose(LDB); 

fclose(CLlCl); 

fclose(CLlC2); 

fclosc(CL2Cl); 

fclose(CL2C2); 

} 

/ * Write a signal of each class to a file for plotting. * / 

void write_signals_to_file(int N) 

{ 

int K; 

double *W0RK1, *W0RK2; 
FILE *P1, *P2; 
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if(( PI = fopcn("/users/stud/eirikf /programs/hovedf ag/idlplot/ 
signalplot/ signall","w")) 
== NULL) { 

priiitf(" ERROR in opening signall\n"); 

exit(l); 

} 

if(( P2 = fopen("/users/stud/eirikf /progreims/hovedf ag/idlplot/ 740 
signalplot/ signal2" ,"w")) 
== NULL) { 

printf("ERROR in opening F2\n"); 

exit(l); 

} 

WORKl = (double *)calloc(N, sizeof (double)); 
W0RK2 = (double *)calloc(N, sizeof (double)); 

niakesignal(WORKl, N, K_l, A_l, CLASS_1); 

makesignal(W0RK2, N, K_2, A_2, CLASS_2); 

fprintf(Pl, "7,d\n", N); 750 
fprintf(P2, "y,d\n", N); 



for (K = 0; K < N; K++) { 
fprintf(Pl, "7.1e\n", W0RK1[K]); 
fprintf(P2, "%le\n", W0RK2[K]); 

} 

frcc(WORKl); 
free(W0RK2); 
fclose(Pl); 
fclose(P2); 

} 



760 



/ * Compute the best decision surface of our two classes of signals 
that can be obtained from the shape of a "hyper-sphere". 

M is the maximum number of coordinates to be used in constructing 
the "hypersphere" . 

1/ 1«P0WER is the step size in the search grid. 

770 

LIMIT is the upper search limit for the best radius. 

LDB is pointer to file containing the coordinates 
of the most discriminating basis elements. 

BESTSPHERE is pointer to file of output data. 

OPTION is a "switch" to reverse the definition of 
the classifier if needed. */ 

780 

void find_decision_surface(int M, int POWER, double LIMIT, 

FILE *LDB, FILE *BESTSPHERE, int OPTION) 
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int N, L, NS, LEVEL, BLOCK, POS, MISPLACED, 
LEASTMISPLACED, BESTDIMENSION, COUNT, LENGTH, J, K, I; 

double *W0RK1, *W0RK2, RADIUS, BESTRADIUS, DELTA, 
STEP, *TRAINP0INT1, *TRAINP0INT2, DISTANCEl, DISTANCE2, 
MISCLASSERROR; 

790 

struct basisnode **LDBDATA; 
FILE *SIG1, *SIG2; 

SIGl = fopen(" /users/ stud/ eirikf /work/timef requeiicy/sigl","r"); 
SIG2 = fopcn("/users/stud/eirikf /work/timef requency/sig2","r"); 
fscanf(LDB, "%d" , &N); 
fscanf(LDB, "7.d", &L); 

fscanfisiGl, "7.d", &N); 

fscanf(SIGl, "7.d", &L); 800 
fscanf(SIGl, "7,d", &NS); 

fscanf(SIG2, , <§^N); 
fscanf(SIG2, "7.d", LL); 
fscanf(SIG2, "7.d", &NS); 
LENGTH = NS*N*(L+1); 

LDBDATA — (struct basisnode **)calloc(M, sizeof (struct basisnode *)); 
WORKl = (double *)calloc(LENGTH, sizeof(double)); 
W0RK2 = (double *)calloc(LENGTH, sizeof(double)); 

for (I = 0; I < LENGTH; I++) { 810 
fscanf(SIGl, "7.1f ", .S^W0RK1[I]); 
fscanf(SIG2, "%lf", &W0RK2[I]); 

} 

TRAINPOINTl = (double *)calloc(M, sizeof(double)); 
TRAINP0INT2 = (double *)calloc(M, sizeof(double)); 

for (K = 0; K < M; K++) { 
fscanf(LDB, "7.d", .S^LEVEL); 
fscanf(LDB, "%d", &BLOCK); 

fscanf(LDB, "7.d", &POS); 820 
fscanf(LDB, "7.1f", &DELTA); 

LDBDATA[K] = makebasisnode(LEVEL, BLOCK, POS, DELTA, 0); 

} 



STEP = 1.0/(1«POWER); 
LEASTMISPLACED = 2*NS; 

for (J = 2; J M; J++) { 

COUNT = 1; 

while (STEP*COUNT <= LIMIT) { 
MISPLACED = 0; 
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RADIUS = STEP*COUNT; 

for (K = 0; K < NS; K++) { 
DISTANCEl = 0; 
DISTANCE2 = 0; 

for (I = 0; I < J; I++) { 
LEVEL = LDBDATA[I]->LEVEL; 

BLOCK = LDBDATA[I]->BLOCK; 840 
POS = LDBDATA[I]->POS; 
TRAINP0INT1[I] = 

*(W0RK1 + K*N*(L+1) + 
abtblock(N,LEVEL,BLOCK) + POS); 
TRAINP0INT2[I] = 

*(W0RK2 + K*N*(L+1) + 
abtblock(N,LEVEL,BLOCK) + POS); 
TRAINPOINTlp] *= TRAINP0INT1[I]; 
TRAINP0INT2[I] TRAINP0INT2[I]; 

DISTANCEl += TRAINP0INT1[I]; 850 
DISTANCE2 += TRAINP0INT2[I]; 

} 

DISTANCEl = sqrt(DISTANCEl); 
DISTANCE2 = sqrt(DISTANCE2); 
if (OPTION 0) { 
if (DISTANCEl < RADIUS) 

MISPLACED += 1; 
if (DISTANCE2 > RADIUS) 
MISPLACED += 1; 
} 860 
else { 

if (DISTANCEl > RADIUS) 

MISPLACED += 1; 
if (DISTANCE2 < RADIUS) 

MISPLACED += 1; 

} 

} 

if (MISPLACED < LEASTMISPLACED) { 

LEASTMISPLACED = MISPLACED; 

BESTDIMENSION = J; 870 
BESTRADIUS = RADIUS; 

} 

COUNT += 1; 

} 

} 

MISCLASSERROR = (1.0*LEASTMISPLACED)/(2.0*NS); 
fprintf(BESTSPHERE, "7.d\ii", BESTDIMENSION); 
fprintf(BESTSPHERE, "7.1f\ii", MISCLASSERROR); 
fprintf(BESTSPHERE, "ZlfXii", BESTRADIUS); 

£ree(WORKl); 880 
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£ree(W0RK2); 

frec(TRAINPOINTl); 
free(TRAINP0INT2); 
freebasisnodes(LDBDATA, M); 
fclose(LDB); 
fclose(BESTSPHERE); 
fclose(SIGl); 
fclose(SIG2); 

} 

/* Given a decision surface to separate our two classes of 
signals, generate a collection of test signals and measure 
the mis classification rate on this collection using 
the given surface. Also compute the misclassification rate 
on the collection of training signals. 

NS is the number of test signals to be generated from each class. 

BESTSPHERE is a pointer to file containing the data 
determining the best classifier using a hyperspherical surface. 900 

LDB is a pointer to file containing the coordinates of 
the most discriminating basis functions. 

DATA is a pointer to the output file. 

OPTION is a "switch" to reverse the definition of the 
classifier. */ 

void test_decision_surface(int NS, FILE *BESTSPHERE, 910 
FILE *LDB, FILE *DATA, int OPTION) 

{ 

int LEVEL, BLOCK, POS, BESTDIMENSION, MISPLACED, 
N, L, J, K, I; 

double *W0RK1, *W0RK2, BESTRADIUS, *TESTP0INT1, 
*TESTP0INT2, MCRTEST, MCRTRAIN, DISTANCEl, 
DISTANCE2, DELTA; 

struct basisnode **LDBDATA; 920 

fscanf (BESTSPHERE, "7.d", &BESTDIMENSION); 
fscanf (BESTSPHERE, "%lf", .S^MCRTRAIN); 
fscanf (BESTSPHERE, "%lf", &BESTRADIUS); 
LDBDATA = (struct basisnode **)calloc(BESTDIMENSION, 

sizeof (struct basisnode)); 

fscanf (LDB, "%d" . .S^N); 
fscanf (LDB, "7.d", &L); 

WORKl = (double *)calloc(N*(L+l), sizeof(double)); 
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W0RK2 = (double *)calloc(N*(L+l), sizeof(double)); 

TESTPOINTl = (double *)calloc(BESTDIMENSION, sizeof (double)); 

TESTP0INT2 = (double *)calloc(BESTDIMENSION, sizeof (double)); 

for (K = 0; K < BESTDIMENSION; K++) { 

fscanf(LDB, "%d" . &LEVEL); 
fscanf(LDB, "%d" , &BLOCK); 
fscanf(LDB, "%d" , &POS); 
fscanf(LDB, "7.1f", &DELTA); 

LDBDATA[K] = makebasisnode(LEVEL, BLOCK, POS, DELTA, 0); 

} 

MISPLACED = 0; 

for (J = 1; J <= NS; J++) { 
niakesignal(WORKl, N, K_l, A_l, CLASS_1); 
makesignal(W0RK2, N, K_2, A_2, CLASS_2); 
expand(WORKl, N, L); 
expand(W0RK2, N, L); 
DISTANCEl = 0; 
DISTANCE2 = 0; 

for (I = 0; I < BESTDIMENSION; I++) { 
LEVEL = LDBDATA[I]->LEVEL; 
BLOCK = LDBDATA[I]->BLOCK; 
POS = LDBDATA[I]->POS; 
TESTPOINTl [I] = 

*(W0RK1 + abtblock(N,LEVEL,BLOCK) + POS); 
TESTP0INT2[I] = 

*(W0RK2 + abtblock(N,LEVEL,BLOCK) + POS); 
TESTPOINTl [I] *= TESTPOINTl [I]; 
TESTP0INT2[I] TESTP0INT2[I]; 
DISTANCEl += TESTPOINTl [I]; 
DISTANCE2 += TESTP0INT2[I]; 

} 

DISTANCEl = sqrt (DISTANCEl); 
DISTANCE2 = sqrt(DISTANCE2); 
if (OPTION == 0) { 
if (DISTANCEl < BESTRADIUS) 

MISPLACED += 1; 
if (DISTANCE2 > BESTRADIUS) 
MISPLACED += 1; 

} 

else { 

if (DISTANCEl > BESTRADIUS) 

MISPLACED += 1; 
if (DISTANCE2 < BESTRADIUS) 
MISPLACED += 1; 

} 

} 
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MCRTEST = (1.0*MISPLACED)/(2.0*NS); 

fpriiitf(DATA, "7.d\n", NS); 

fprintf(DATA, "7.d\n", BESTDIMENSION); 

fprintfpATA, "7.1f\n", MCRTRAIN); 

fprintf(DATA, "7,lf\n", MCRTEST); 

fprintfpATA, "7,lf\n", BESTRADIUS); 

free(WORKl); 

£ree(W0RK2); 

free(TESTPOINTl); 

frcc(TESTP0INT2); 

freebasisnodes(LDBDATA, BESTDIMENSION); 

fclose(LDB); 

fclose(BESTSPHERE); 

fclose(DATA); 
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3. Author's Maple code 



# Compute the S02(R) elements from the filter coefficients. # 

lagevinkler := 
proc (k) ; 

for i from 1 to N do 

if i < N then x(i) := evalf ( h(i-l,i-l) / h(i-l,i) ) fi ; 
for j from i to 2 * N — i do 
if type ( i + j, odd ) then 
h(i,j) := evalf ( - x(i) * h(i-l,j+l) + h(i-l,j) ) 
else 10 
h(i,j) := evalf ( x(i) * h(i-l,j-l) + h(i-l,j) ) 
fi ; 

od ; 

if i = N then x(N) := evalf ( - h(N-l,N) / h(N-l,N-l) ) 

fi ; 
od ; 

for i from 1 to N do print ( i, x(i) ) 
od ; 
end ; 

20 

# Compute the polynomials P{n}(m) of degree n that result from lowpass 
filtering of the monomials m"{n}. # 



f := proc( j,i,m ) ; 

if m = then RETURN(1/10) else 
RETURN( ( (j+2*i)/10)-m ) fi ; 
end ; 

30 



z(l) 


= .4122865951 




z(2) 


1.831178514 ; 


z(3) 


= -.1058894200 ; 


z(4) 


= 7.508378888 




z(5) 


-.02083494630 


z(6) 


= -.04396989341 


z(7) 


= -.004543409641 


z(8) 


= 15.03636438 




z(9) 


= -.009120773147 


z(10) 


:= 4.100416655 




z(ll) 


:= 15.61099604 




z(12) 


:= 11.59905847 




z(13) 


:= 37.56973541 




z(14) 


:= 197.1316159 




z(15) 


:= -.0004543371650 
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N — 3 ; 

50 

javel := proc( 1 ) 

for m from to 6 do 
for i from —4 to 4 do 

for j from to 2*N— 1 do a(m,0,j,i) := evalf( f(j,i,m) ) od ; 

for k from 1 to N do 

for j from k— 1 to 2*N— k do 

if type ( k+j,odd ) then 

a(m,k,j,i) :— evalf( — z(k)*a(m,k— l,j+l,i)+a(m,k— l,j,i) ) else 60 
a(m,k,j,i) := evalf( z(k)*a(m,k— l,i)+a(m,k— l,j,i) ) 

fi ; 
od ; 
od ; 
od ; 
od ; 



pO := a(0,N,N-l,0) ; 
pi := n -> Bl*n + CI ; 
p2 := n -> B2*n"2 + C2*n + D2 ; 
p3 
p4 

s := solvc( { pl(N-l)=a(l,N,N-l,0), 
pl(N+l)=a(l,N,N-l,l) },{ Bl.Cl } ) ; 
assign( s ) ; 

s := solve( { p2(N-l)=a(2,N,N-l,0), 
p2(N+l)=a(2,N,N-l,l), 
p2(N+3)=a(2,N,N-l,2) },{ B2,C2,D2 } ) ; 
assign ( s ) ; 
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:= n -> B3*n"3 + C3*n-2 +D3*n + E3 ; 
:= n -> B4*n"4 + C4*n"3 + D4*n"2 + E4*n + F4 ; 



s := solve( { p3(N-3)=a(3,N,N-l,-l), 

p3(N-l)=a(3,N,N-l,0),p3(N+l)=a(3,N,N-l,l), 

p3(N+3)=a(3,N,N-l,2) },{ B3,C3,D3,E3 } ) ; 90 
assign ( s ) ; 

s := solve( { p4(N-5)=a(4,N,N-l,-2), 
p4(N-3)=a(4,N,N-l,-l),p4(N-l)=a(4,N,N-l,0), 

p4(N+l)=a(4,N,N-l,l),p4(N+3)=a(4,N,N-l,2) },{ B4,C4,D4,E4,F4 } ) ; 
assign( s ) ; 
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print ( pO ) ; 
end ; 

# Compute the row vectors of the edge matrices by claiming 
some vanishing moments on the highpass edge coefficients and 
certain polynomial conditions on the lowpass edge coefficients ensuring 
polynomials of degree up to n mapping continuously to polynomials of 
the same degree. # 

with( linalg ) ; 

jafs := proc( o ) 

v(0) array( 1. .5, [a(0, 0,0,0), a(0, 1,1,0), a(0, 2,2,0), a(0, 3,3,0), 0] ) ; 
v(l) := array( 1. .5,[a(l,0,0,0),a(l,l,l,0),a(l,2,2,0),a(l,3,3,0),0] ) ; 
v(2) := array( 1. .5,[a(2,0,0,0),a(2,l,l,0),a(2,2,2,0),a(2,3,3,0),0] ) ; 
v(3) := array( 1. .5, [a(3, 0,0,0), a(3, 1,1,0), a(3, 2,2,0), a(3, 3,3, 0),0] ) ; 
v(4) := array( 1. .5,[a(4,l,0,0),a(4,2,l,0),a(4,3,2,0),a(4,4,3,0), 
a(4,5,4,0)] ) ; 



cO :— array( 1. .5,[] ) ; 
c0[5] := ; 

do := array( 1. .5,[] ) ; 
do [4] :=0 ; 
d0[5] := ; 

cl := array( 1. .5,[] ) ; 

dl := array( 1. .5,[] ) ; 
dl[5] := ; 

d2 := array( 1. .5,[] ) ; 

s := solvc( { dotprod( cO,v(0) ) = pO, dotprod( c0,v(l) ) = pl(l), 
dotprod( c0,v(2) ) = p2(l), dotprod( c0,v(3) ) = p3(l) }, 
{ c0[l],c0[2],c0[3],c0[4] } ) ; 
assign( s ) ; 

s := solve( { dotprod( dO,v(l) ) = 0, dotprod( dO,v(0) ) = }, 
{ d0[l].d0[2] } ) ; 
assign( s ) ; 

s := solve( { dotprod( cl,v(0) ) = pO, dotprod( cl,v(l) ) = pl(3), 
dotprod( cl,v(2) ) = p2(3), dotprod( cl,v(3) ) = p3(3) }, 
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{ cl[l],cl[2],cl[3],cl[4] } ) ; 
assign( s ) ; 

s := solve( { dotprod( dl,v(l) ) = 0, dotprod( dl,v(2) ) = 0, 
dotprod( dl,v(0) ) = }, 
{ dl[l].dl[2],dl[3] } ) ; 
assign( s ) ; 



iO := -10.5 ; 
jO := -10.5 ; 
10 := -10.5 ; 

d0[3] := iO ; 
cl[5] := jO ; 
dl[4] := 10 ; 

A evalm( array( 1. .5,1. .5, [(1,1)= dG[l],(l,2)=dO[2],(l,3)=dG[3], 

(1.4) =d0[4],(l,5)=d0[5],(2,l)=c0[l],(2,2)=c0[2],(2,3)=c0[3],(2,4)=c0[4], 

(2.5) =c0[5],(3,l)=dl[l],(3,2)=dl[2],(3,3)=dl[3],(3,4)=dl[4],(3,5)=dl[5], 

(4.1) =cl[l],(4,2)=cl[2],(4,3)=cl[3],(4,4)=cl[4],(4,5)=cl[5],(5,l)=0, 

(5.2) =0,(5,3)=0,(5,4)=0,(5,5)=1 ] ) ) ; 

B := evalni( inverse( A ) ) ; 

z := evalf( norni( B,l ) + norm( A,l ) ) ; 

for i from iO by 1 to 10.5 do 
for j from jO by 1 to 10.5 do 
for 1 from jO by 1 to 10.5 do 

d0[3] := i ; 
cl[5] j ; 
dl[4] := 1 ; 

A := evalm( array( 1. .5,1. .5, [(1,1)= d0[l],(l,2)=d0[2],(l,3)=d0[3], 

(1.4) =d0[4],(l,5)=d0[5],(2,l)=c0[l],(2,2)=c0[2],(2,3)=c0[3],(2,4)=c0[4], 

(2.5) =c0[5],(3,l)=dl[l],(3,2)=dl[2],(3,3)=dl[3],(3,4)=dl[4],(3,5)=dl[5], 

(4.1) =cl[l],(4,2)=cl[2],(4,3)=cl[3],(4,4)=cl[4],(4,5)=cl[5],(5,l)=0, 

(5.2) =0,(5,3)=0,(5,4)=0,(5,5)=1 ] ) ) ; 

B := evalm( inverse( A ) ) ; 

w := evalf( norm( B,l ) + norm( A,l ) ) ; 



if w < z then z 



= w 
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fi ; 

od 
od 
od 



d0[3] := iO ; 
cl[5] := jO ; 
dl[4] := 10 ; 

A cvalm( array( 1..5,1..5, [(1,1)= d0[l],(l,2)=d0[2],(l,3)=d0[3], 

(1.4) =d0[4],(l,5)=d0[5],(2,l)=c0[l],(2,2)=c0[2],(2,3)=c0[3],(2,4)=c0[4], 

(2.5) =c0[5],(3,l)=dl[l],(3,2)=dl[2],(3,3)=dl[3],(3,4)=dl[4],(3,5)=dl[5], 

(4.1) =cl[l],(4,2)=cl[2],(4,3)=cl[3],(4,4)=cl[4],(4,5)=cl[5],(5,l)=0, 

(5.2) =0,(5,3)=0,(5,4)=0,(5,5)=1 ] ) ) ; 

B := evalni( inverse( A ) ) ; 
print ( z ) ; 

print( l,norni( A,l ) ) ; 
print( 2,norni( B,l ) ) ; 
print( 1,A ) ; 
print( 2,B ) ; 

print ( evalm( B A ) ) ; 



end 
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